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Once  a  uniform  electric  field  is  turned  on  in  graphene,  carriers  accelerate  ballistically  until  they  are  scattered 
by  optic  phonons  and  the  process  repeats  itself.  In  this  dissertation,  I  will  show  that  the  oscillatory  nature 
of  the  motion  of  the  carrier  distribution  function  manifests  in  damped  oscillations  of  carrier  drift  velocity 
and  average  energy.  In  appropriate  fields,  the  frequency  of  such  oscillations  can  be  in  the  terahertz  (THz) 
range.  The  randomizing  nature  of  optical  phonon  scattering  on  graphene’s  linear  band  structure  further 
limits  terahertz  observation  to  a  range  of  sample  lengths. 

I  will  also  show  that  when  an  ac  field  is  superimposed  onto  the  appropriate  dc  field,  hot  carriers  in 
graphene  undergo  an  anomalous  parametric  resonance.  Such  resonance  occurs  at  about  half  the  frequency 
?F  =  2?cF/~?OP  ,  where  2?/?F  is  the  time  taken  for  carriers  to  accelerate  ballistically  to  the  optic  phonon 
energy  ~?OP  in  a  dc  field  F.  For  weak  elastic  scattering,  the  phase  difference  between  the  current  and  the 
ac  field  has  a  nonzero  minimum  at  resonance.  Dephasing  increases  with  ac  frequency  for  stronger  elastic 
scattering.  The  overall  effect  could  also  be  seen  in  long-range  spatially  periodic  potentials  under  steady  state 
conditions. 

This  dissertation  also  shows  that  the  soft  parametric  resonance  (SPR)  at  ?  =  ??F  is  temperature 
independent,  and  the  resonance  factor  ?  ?  0.56  is  weakly  dependent  on  the  dc  field  Fo.  This  ensures 
tunability  of  resonant  frequencies  in  the  terahertz  range  by  varying  Fo.  A  small  signal  analysis  (SSA)  of  the 
time-dependent  Boltzmann  transport  equation  (BTE)  reveals  a  second  resonance  peak  at  ?  ?  1 .  This  peak 
is  prevalent  at  temperatures  T  ?  77  K,  and  appears  as  a  weak  shoulder  at  T  =  300  K. 

Finally,  this  dissertation  shows  that  in  graphene,  the  motion  of  carriers  under  the  influence  of  temporarily 
and  spatially  modulated  scattering  is  characterized  by  sharp  resonances.  Such  resonances  occur  when  the 
period  of  the  ac  field  applied  equals  the  time  taken  by  the  quasi-ballistic  carriers  to  travel  a  spatial  distance 
corresponding  to  the  wavelength  of  the  field.  I  will  also  show  that  such  scattering  can  be  realized  on  graphene 
sheets  on  periodically  spaced  gates  energized  by  an  a-c  bias.  Appropriate  fields  and  gate  separation  will 
result  in  high  Q-factor  resonances  in  the  THz  range.  The  resonant  frequencies  are  tunable  with  the  gate 
separation,  and  higher  harmonics  with  large  Q-factors  can  also  be  achieved. 
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Abstract 


Once  a  uniform  electric  field  is  turned  on  in  graphene,  carriers  accelerate  ballistically  until  they  are  scattered 
by  optic  phonons  and  the  process  repeats  itself.  In  this  dissertation,  I  will  show  that  the  oscillatory  nature 
of  the  motion  of  the  carrier  distribution  function  manifests  in  damped  oscillations  of  carrier  drift  velocity 
and  average  energy.  In  appropriate  fields,  the  frequency  of  such  oscillations  can  be  in  the  terahertz  (THz) 
range.  The  randomizing  nature  of  optical  phonon  scattering  on  graphene’s  linear  band  structure  further 
limits  terahertz  observation  to  a  range  of  sample  lengths. 

I  will  also  show  that  when  an  ac  field  is  superimposed  onto  the  appropriate  dc  field,  hot  carriers  in 
graphene  undergo  an  anomalous  parametric  resonance.  Such  resonance  occurs  at  about  half  the  frequency 
ujf  =  2neF/hu)oPi  where  2 n/uip  is  the  time  taken  for  carriers  to  accelerate  ballistically  to  the  optic  phonon 
energy  Hu>op  in  a  dc  field  F.  For  weak  elastic  scattering,  the  phase  difference  between  the  current  and  the 
ac  field  has  a  nonzero  minimum  at  resonance.  Dephasing  increases  with  ac  frequency  for  stronger  elastic 
scattering.  The  overall  effect  could  also  be  seen  in  long-range  spatially  periodic  potentials  under  steady  state 
conditions. 

This  dissertation  also  shows  that  the  soft  parametric  resonance  (SPR)  at  to  =  r)LUp  is  temperature 
independent,  and  the  resonance  factor  77  ~  0.56  is  weakly  dependent  on  the  dc  field  Fa.  This  ensures 
tunability  of  resonant  frequencies  in  the  terahertz  range  by  varying  F0.  A  small  signal  analysis  (SSA)  of  the 
time-dependent  Boltzmann  transport  equation  (BTE)  reveals  a  second  resonance  peak  at  77  ~  1.  This  peak 
is  prevalent  at  temperatures  T  <  77  K,  and  appears  as  a  weak  shoulder  at  T  =  300  K. 

Finally,  this  dissertation  shows  that  in  graphene,  the  motion  of  carriers  under  the  influence  of  temporarily 
and  spatially  modulated  scattering  is  characterized  by  sharp  resonances.  Such  resonances  occur  when  the 
period  of  the  ac  field  applied  equals  the  time  taken  by  the  quasi-ballistic  carriers  to  travel  a  spatial  distance 
corresponding  to  the  wavelength  of  the  field.  I  will  also  show  that  such  scattering  can  be  realized  on  graphene 
sheets  on  periodically  spaced  gates  energized  by  an  a-c  bias.  Appropriate  fields  and  gate  separation  will 
result  in  high  Q-factor  resonances  in  the  THz  range.  The  resonant  frequencies  are  tunable  with  the  gate 
separation,  and  higher  harmonics  with  large  Q-factors  can  also  be  achieved. 
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Chapter  1 


INTRODUCTION 


1.1  TERAHERTZ  RADIATION 

Terahertz  (THz)  radiation  is  electromagnetic  radiation  with  frequencies  in  the  range  0.3-3  THz.1  This  is 
a  region  in  the  electromagnetic  spectrum  between  the  microwave  and  infrared  regions.  Because  THz  waves 
have  wavelengths  below  1  mm,  they  are  also  referred  to  as  submillimeter  waves.  This  remains  the  least 
developed  region  in  the  spectrum  because  of  the  lack  of  convenient  sources  and  detectors  of  THz  waves. 
Consequently,  this  region  of  the  spectrum  is  also  known  as  the  “terahertz  gap” . 2  THz  waves  have  attracted 
a  lot  of  interest  in  the  last  two  decades  because  of  their  wide  range  of  potential  applications.3,4,5,6,7,8  THz 
radiation  could  be  useful  in  astrophysics  and  atmospheric  science,  biological  and  medical  science,  security 
screening,  non-destructive  evaluation,  communications  technology,  and  ultrafast  spectroscopy. 9,10,11 

Terahertz  radiation  actually  occurs  naturally  as  black-body  radiation  from  objects  with  temperatures 
grater  than  10  K.12  Even  though  this  thermal  emission  is  very  weak  for  most  practical  uses,  it  is  still 
useful  in  astrophysics. 13  Terahertz  waves  can  also  be  generated  using  transistor  oscillators  and  amplifiers, 
but  the  output  power  is  too  low  for  most  applications.  The  power  generated  by  such  devices  falls  off 
at  frequencies  higher  than  microwaves  as  a  result  of  transit-time  and  resistance-capacitance  effects. 14,15 
Designing  conventional  semiconductor  diode  lasers-  that  generate  infrared  and  visible  radiation-to  generate 
THz  radiation  has  also  proved  difficult.  This  is  because  of  the  lack  of  materials  with  small  enough  band 
gaps  to  generate  electromagnetic  waves  with  frequencies  below  15  THz.9 

Significant  research  efforts  has  led  to  the  development  of  some  solid-state  THz  devices.4  Other  sources 
of  THz  radiation  include  THz  quantum  cascade  lasers  and  resonant  tunneling  diodes.9,16,17,18,19  Despite 
significant  improvements  in  the  performance  of  such  devices,20,1,21  there  is  still  a  need  for  coherent,  compact, 
high-power,  and  tunable  THz  sources  and  detectors  that  operate  at  room  temperature. 
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1.2  GRAPHENE 


Graphene  is  a  two-dimensional  (2D)  material  with  unique  electronic  properties. 22  Even  though  the  material 
had  been  extensively  studied  many  years  before,23  it  was  not  until  2004  that  graphene  was  isolated.24  This 
was  a  remarkable  turn  of  events  because  it  was  widely  known  that  strictly  2D  materials  could  not  exist.25,26 
The  isolation  of  graphene  led  to  “graphene  rush”  as  the  material  became  the  focus  of  extensive  fundamental 
and  technological  studies  by  condensed  matter  physicists,  chemists,  and  engineers.  What  really  sets  graphene 
apart  from  conventional  semiconductor  is  that  it’s  electrons  behave  as  massless  charged  fermions  and  can  be 
described  by  a  formalism  similar  to  the  Dirac  relativistic  equation  rather  than  the  Schroendinger  equation. 2 ' 
This  presents  an  opportunity  to  study  quantum  electrodynamics  at  the  nanoscale  level.28,29  Other  exciting 
properties  of  graphene  include  excellent  electrical  conductivity,  optical  transparency,  mechanical  strength, 
and  chemical  stability. 30 


Figure  1.1:  (Color  online)  Left:  Graphene  lattice  with  two  atoms  per  unit  cell.  The  vectors  ai  and  a2  are  the 
unit  vectors  and  nearest-neighbor  vectors  are  given  by  5i,i  =  1,2,3.  Right:  the  corresponding  Brillouin  zone.  The 
conduction  and  valence  bands  in  graphene  meet  at  the  points  K  and  K'  called  Dirac  points.  From  Castro  Neto, 
Guinea,  Peres,  Novoselov,  and  Geim,  Rev.  Mod.  Phys.  81,  109  (2009).  27 


1.2.1  Graphene  Band  Structure 

Graphene  is  made  up  of  carbon  atoms  arranged  in  a  honeycomb  structure  (Fig.  1.1).  Each  unit  cell  of 
graphene  consists  of  two  carbon  atoms  and  is  repeated  along  the  lattice  vectors  R  =  mdi  +  nai) ,  where 

a!  =av/3(i,-^),  a2  =av/3(^,-^)  (1.1) 
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and  a  ~  1.42 A  is  the  inter-atomic  distance.  The  corresponding  Brillouin  zone  (BZ)  has  the  lattice  vectors 
(Fig.  1.1) 

-*  9tT  -*  O'TT 

bl  =  —(l,V3),  b2  =  — (1,  —VS).  (1.2) 

da  da 

The  corner  points  K  and  K'  of  the  BZ  play  an  important  role  in  the  transport  properties  of  graphene.  Their 
locations  in  momentum  space  given  by 


T-j .  2tt  1  2n ,  1 

3^(  =  3^( 


A  tight-binding  model  gives  the  Hamiltonian  of  the  electrons  in  graphene  as 


(1.3) 


H  =  -ht  ^2  +  h.C,  (1.4) 

ra,n,er 

where  m(aa  m)  creates  (annihilates)  an  electron  with  spin  a(a  =  (+,—))  on  site  Rm  on  sublattice  A,  and 
h.c  is  the  hermitian  conjugate  of  the  first  term.  The  parameter  t  ss  2.8  eV  is  the  hopping  energy  between 
two  nearest  lattice  sites.  In  the  limit  of  low  energies  (more  relevant  for  electronic  transport),  hoppings  to 
next  nearest  and  further  neighbours  do  not  play  a  role  and  have  not  been  included  in  Eq.  1.4.  The  energy 
bands  corresponding  to  the  Hamiltonian  (1.4)  are  given  by23 


E(k)  =  ±Hty  3  +  /(fc), 

f(k)  =  2  cos  (VSkya)  +  4  cos  (-^- kya )  cos  (- kxa ), 


(1.5) 


where  the  plus  sign  is  for  the  upper  (n*)  and  the  minus  sign  is  for  the  lower  (7r)  band  (Fig.  1.2).  Note 
that  the  two  bands  meet  at  six  zero  energy  points  called  Dirac  or  neutrality  points.  By  symmetry,  these 
points  correspond  to  the  two  independent  points  K  and  K'  in  the  BZ  of  graphene.22,31  Because  the  two 
bands  touch  at  the  Dirac  points,  graphene  has  no  band  gap,  and  is  referred  to  as  a  zero-gap  semiconductor. 
Expanding  Eq.  1.5  around  K  or  ( I\ ')  we  get,23 


E(k)  ~  ±Hvp\k\  (1.6) 

for  k  close  to  K(K'),  and  Vf  =  3fa/2  «  1  x  108  cm/s.  As  a  result,  carriers  in  graphene  travel  with  Fermi 
velocity  vf  that  is  much  larger  than  in  conventional  semiconductors. 32  The  high  value  of  vf  and  superior 
conductivity  make  graphene  an  attractive  material  for  future  high  speed  electronic  devices.  33.34>35>36>37 
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Figure  1.2:  (Color  online)  Graphene  band  structure.  Inset:  Band  structure  close  to  one  of  the  Dirac  points.  From 
Castro  Neto,  Guinea,  Peres,  Novoselov,  and  Geim,  Rev.  Mod.  Phys.  81,  109  (2009).  27 

Because  there  are  two  sublattices  A  and  B  in  the  structure  of  graphene  (Fig.  1.1),  the  Hamiltonian 
approximates  to  the  Dirac  equation 

H  =  hvFcr  ■  k,  (1.7) 

where  er  represents  pseudospin  due  to  the  two  atoms  in  the  unit  cell.  Thus,  electrons  in  graphene  behave 
like  charged  relativistic  particles  that  have  zero  rest- mass  and  move  with  speed  vf-  As  a  result,  graphene 
provides  a  platform  to  study  the  physics  of  massless  Dirac  fermions  in  solids. 

1.3  OVERVIEW 

The  main  goal  of  this  dissertation  is  to  highlight  some  of  the  properties  of  THz  oscillations  of  hot  electrons 
in  graphene.  Chapter  2  contains  a  discussion  of  dynamics  of  electrons  in  the  conduction  band  of  graphene, 
when  a  a  uniform  electric  field  is  applied.  The  chapter  also  discusses  the  necessary  conditions  required  to 
produce  room  temperature  THz  oscillations  of  the  electron  distribution  in  graphene. 

Chapter  3  includes  a  study  of  parametric  oscillations  of  hot  carriers  in  graphene  under  the  influence  of 
an  ac  field  superimposed  onto  a  dc  field.  Also  discussed  are  the  conditions  for  occurrence  of  parametric 
resonance  in  the  THz  range,  and  the  tunability  of  such  frequencies  with  the  dc  field.  Chapter  4  discusses 
the  effects  of  temperature  and  the  dc  field  on  the  THz  parametric  resonance  condition.  A  discussion  on  a 
secondary  resonance  peak  that  occurs  at  low  temperatures  is  also  included  in  the  chapter. 

The  study  of  electron  dynamics  in  graphene  under  the  influence  of  spatially  and  temporary  modulated 
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electric  fields  is  included  in  chapter  5.  Such  electric  fields  and  scattering  can  be  realized  by  placing  free 
standing  graphene  sheets  on  periodically  spaced  gates  modulated  by  an  ac  field.  The  chapter  discusses 
the  nature  of  the  high  Q-factor  THz  resonances  observed,  the  generation  of  higher  harmonics,  and  their 
t  unability. 

The  last  chapter  discusses  the  important  results  of  my  work  and  the  implications  in  THz  science  and 
technology. 
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Chapter  2 


HOT-ELECTRON  TRANSIENT 
AND  TERAHERTZ  OSCILLATIONS 
IN  GRAPHENE* 

2.1  INTRODUCTION 

In  graphene,  the  carrier  mean  free  path  for  collisions  with  low  energy  acoustic  phonons  (AP)  can  reach  several 
micrometers,  while  efficient  scattering  with  monochromatic  optic  phonons  (OP)  occurs  at  much  larger  energy 
( tvujop  ~  0.2  eV).38’39  Therefore,  in  electric  fields  F  high  enough  (to  escape  AP  or  low-energy  scattering), 
charge  carriers  can  experience  quasiballistic  runaway  until  they  scattered  with  efficient  OP  emission  once 
E  >  Huj0p-40  The  coherent  aspect  of  this  process,  that  is,  the  coherent  acceleration  of  the  carrier  distribution 
function  followed  by  quasi-instantenous  carrier  relaxation  by  OP’s  at  high  energy,  is  expected  to  produce 
oscillations  of  the  carrier  velocity  with  a  periodicity  given  by  Tgr  =  huop / eFvF  ~  1  ps  (Ref.  41)  for  F  ~  2 
kV/cm,  that  is,  in  the  terahertz  range. 

This  type  of  oscillation  has  already  been  predicted  in  GaAs,41  and  indirectly  observed  in  slightly  n- 
doped  InSb.42  However,  their  manifestation  in  these  materials  is  different  in  several  respects:  First,  owing 
to  the  carrier  parabolic  energy-momentum  dispersion,  the  oscillation  periodicity  in  GaAs  is  instead  given  by 
t GaAs  =  \J 2to  *  HloqpAs/ eF .  41,43,44  geconcj;  and  more  importantly,  in  III-V  semiconductors,  the  oscillation 
onset  is  restricted  by  two  conflicting  conditions:  On  the  one  hand,  the  low  value  of  the  OP  energy  ( Hujop  ~ 
0.04  eV)  (Ref.  45)  is  comparable  to  the  thermal  broadening  of  the  carrier  distribution  at  room  temperature 
so  that  the  back-and-forth  motion  of  the  distribution  between  the  optic  phonon  and  the  zero  point  energy 
is  immediately  damped. 41,43,44  On  the  other  hand,  at  low  temperature,  ionized  impurity  scattering  becomes 
dominant  and  produces  strong  damping  which  can  only  be  reduced  by  lowering  the  dopant  density,  thereby 
lowering  the  carrier  density,  and  weakening  the  oscillation  amplitude.  In  this  respect,  the  high  conductance 
of  graphene,  and  the  high  optic  phonon  energy  provide  the  conditions  for  room-temperature  observation. 

This  chapter  provides  a  study  of  carrier  dynamics  in  graphene  single  layers,  once  the  electric  field  has 
been  turned  on,  to  determine  the  conditions  of  occurrence  of  current  oscillations.  Indeed,  thermal  broadening 
of  the  initial  distribution  introduces  a  decoherence  in  the  OP  relaxation  among  carriers,  which,  aside  from 

*This  chapter  closely  follows  Samwel  Sekwao,  and  Jean-Pierre  Leburton,  Phys.  Rev.  B  83,  075418  (2011). 
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low-energy  scattering,  causes  inherent  damping  of  the  oscillations,  especially  at  high  fields,  where  carriers 
overshoot  the  OP  energy  before  emission.  For  this  purpose,  the  Boltzmann  formalism  is  used  to  solve 
for  the  steady  state  and  the  time-dependent  carrier  distribution  in  the  presence  of  OP  scattering.  41,43>44 
Since  both  intravalley  and  intervalley  scattering  can  be  treated  within  the  deformation-potential  interaction 
model,46  an  overall  transition  rate  accounting  for  both  mechanisms  is  used.  Low-energy  scattering  such 
as  impurities  and  acoustic  phonons  by  are  accounted  for  using  the  relaxation  time  approximation, 47  but 
the  dissipation  effects  due  to  the  environment  such  as  remote  phonons,48  detrimental  to  the  occurrence  of 
the  oscillations  are  ignored,  which  limits  their  observations  optimally  to  suspended  graphene  or  nonpolar 
substrate.49  In  weak  concentrations  (nc  ~  1011  cm-2),  electron-electron  interactions  do  not  play  a  major 
role  in  transport  in  graphene,27  and  they  are  not  included  in  this  analysis.  This  chapter  also  analyses  the 
interplay  between  applied  electric  fields  and  strength  of  the  low-  energy  scattering  rates  for  the  onset  of 
oscillations  is  provided.  In  particular  it  is  shown  that  unlike  in  III-V  semiconductors,  where  the  OP  polar 
nature  focuses  the  low-energy  repopulation  along  the  field,  and  provides  a  “streaming”  profile  to  the  carrier 
distribution,  in  graphene,  the  randomizing  scattering  by  deformation  potential  OP  generates  of  a  transient 
crater-like  shape  in  the  low-energy  distribution. 


2.2  MODEL 


Let  us  consider  a  system  of  electrons  in  the  graphene  conduction  band  under  the  influence  of  an  electric  field 
along  the  x-direction.  The  momentum  space  can  be  divided  into  two  regions:  I  (fc  <  kc )  and  II  ( k  >  kc ) 
separated  by  a  circle  of  critical  momentum  kc  =  ujop/vf ,  which  corresponds  to  the  electron  kinetic  energy 
equal  to  fnoop  [Fig.  2.1(a)].  In  region  I  electrons  undergo  quasiballistic  acceleration  and  weak  scattering  by 
low-energy  mechanisms  (e.g.,  impurities  or  AP’s)  until  they  reach  region  II  where  they  lose  their  energy  by 
OP  emission,  and  scatter  back  to  region  I.  In  this  model,  the  electric  field  is  assumed  to  be  low  enough  so 
that  electrons  are  scattered  efficiently  from  region  II  to  region  I  by  OP  emission  with  little  probability  to 
reach  E  >  2 Hojop- 

Quite  generally,  the  time-dependent  Boltzmann  equation  in  each  region  can  be  written  as; 


dfijk,  t)  eF  df!(k,t ) 

dt  H  dkx 


h(k,t)-f0(k) 

r 


(2.1a) 
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Figure  2.1:  (Color  online)  a)  Schematics  of  carrier  quasiballistic  acceleration  and  OP  scattering  in  2D  k  space.  The 
solid  circle  of  radius  kc  is  the  locus  of  all  points  in  k  space  corresponding  to  the  carrier  energy  huiop-  b)  Schematics 
of  electrons  scattered  by  OP’s  from  two  different  positions  of  the  distribution  function  shown  by  ovals  in  region  II. 
The  dashed  circles  in  region  I  represent  the  two  areas  of  high  momentum  probability  where  electrons  are  more  likely 
to  land  after  OP  emission  from  the  distribution  in  region  II  at  two  different  times.  The  arrows  represent  the  electron 
transitions  from  II  to  I. 


(2.1b) 


dfu(k,  t )  eF  dfu(k,t) 

dt  h  dkx 


-fn(k,t)Y,S(k,k') 

k' 


where  fi(k,t)  and  fu(k,t)  are  the  time-dependent  momentum  distribution  functions  in  regions  I  and  II, 
respectively,  and  F  is  the  electric  field  applied.44  The  first  term  on  the  right-hand-side  (RHS)  of  Eq.  (2.1a) 
accounts  for  low-energy  scattering  mechanisms  where, 


fo(k)  = 


l  +  exp(^-fcF)) 


(2.2) 


is  the  Fermi-Dirac  equilibrium  distribution  function  and  fcp  >  0  is  the  Fermi  momentum  so  that  only 
intraband  processes  restricted  to  the  conduction  band  are  considered.  Interband  processes  will  be  the  object 
of  future  work.  Also,  note  that  at  room  temperature,  for  Fermi  levels  at  Ep  =  kpT  above  the  Dirac  point, 
the  carrier  concentrations  is  nc  ss  1.8  x  1011  cm'2.  The  parameter  r  is  the  relaxation  time50  (for  the  sake 
of  simplicity,  it  is  assumed  that  r  is  /c-independent  and  it’s  value  is  varied  compared  to  the  OP  scattering 
rate).  S(k',k)  is  the  OP  transition  rate  from  a  state  with  momentum  k'  to  the  state  with  momentum  k, 
given  by f  51,52 

_  _  ^ n  2n\r  b1  _l 

(2.3) 


an  rn  ttD0  (Nq  +  1)8(E'  -E+  Hujop) 


where  Da  is  an  effective  optical  deformation  potential  accounting  for  intra  and  intervalley  scattering.  The 
parameter  er  is  the  mass  of  the  graphene  sheet  per  unit  area,  and  A  is  the  area  of  the  sheet. 

The  second  term  on  the  (RHS)  of  Eq.  (2.1b)  is  the  carrier  depopulation  by  OP  emission,  while  in 
Eq.  (2.1a)  it  is  the  corresponding  carrier  repopulation  at  low  energy.  Here  OP  absorption  processes  are 
neglected  since  for  temperature  T  =  300  K  used  throughout  this  analysis,  fnoop  ^  kpT,  for  which  the 
phonon  occupation  number  Nq  is  negligible,  so  OP’s  only  scatter  electrons  from  region  II  to  region  I. 


2.3  STEADY  STATE  REGIME 


First,  the  solution  for  //  and  ///  under  steady-state  conditions  is  obtained  by  setting  =  0 

in  Eqs.  (2.1a)  and  (2.1b).  The  procedure  is  to  solve  Eq.  (2.1b)  for  ///,  and  substitute  it’s  value  in  Eq. 
(2.1a)  to  solve  for  //.  Then  one  matches  the  two  solutions  at  the  boundary  k  =  kc ■ 

f  The  transition  rate  should  be  multiplied  by  (1  +  008(0')),  where  0'  is  the  angle  between  k  and  k' .  The  summation  over 
cos(0')  prohibits  scattering  to  areas  near  —kx  axis  in  region  I.  This,  however  does  not  introduce  significant  changes  to  the 
results. 
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The  general  solution  of  Eq.  (2.1b)  is  then, 


fii{kx,kv)  =  fi{k°x,ky)T{kx,ky) 


where  fi(kx,ky)  =  fn(kx,ky)  is  the  distribution  function  at  the  boundary  and  kx  =  ■ \J kc2  —  ky2 
value  at  the  boundary  k  =  kc.  From  here  on,  the  transformation  kx/kc  — >  kx,  ky/kc  — >  ky ,  and  kc 
be  used.  In  this  framework,  the  function  T  is  given  by, 


T(fcx,  ky)  =  exp  <  -a  /  {\I z2  +  kv  -  l}dz 


lk° 


where  the  dimensionless  parameter  a  is  given  by 


Po2kc2 

2auJopeFvF 


Substituting  Eq.  (2.4)  into  Eq.  (2.1a)  leads  to  the  following  equation  in  region  I: 


9fi  ,  *  tn  ,  \  ,  a(k+l)  f1  T(kxc,y) 

+  aiJi  =  “7 Jo{kx,  kv)  +  — — —  J  i  dyfbiv) — ^ - , 


where 


kxc  =  V(k+  l)2  -  y 2 


with 


k  =  \!  kx  +  ky 


and  the  damping  parameter  7  is  defined  as 


To 

7  =  —  > 


where  rQ  =  rop(k  =1.5 kc)  and 


rop(k) 


=  E^^') 


is  the  OP  scattering  rate.  Using  tyjjop  =  0.2  eV,  Da  =  14  eV/A,53  and  <r  =  7.61  x  10  7  kg/m2, 
t0  Ri  0.32ps.  The  preceding  equation  is  a  first-order  partial  differential  equation  with  solution 


CL  f  ^ 

f  1  ( kx ,  ky)  1 1  ( kx ,  ky)  +  -  I  dyfb(y)K  {kx ,  y ,  ky ) , 


(2.4) 

is  the  kx 

1  will 


(2.5) 


one  finds 

(2.6) 
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where 


and 


with 


Tx{kx,  ky)  =  a7e-Q^  J  ^  dzfD(z:  ky)e^z 

fkx  dzT(kxc{z),y){Jz2  +  k2  +  l}e“7Z 

K(kx,  y,  ky)  =  e~ajkx  /  - -V— - 

Jk0  kxc{z) 

kxc{z)  =  \J{yJz2  +  k2  +  l)2-y2. 


Evaluating  fj  at  the  boundary  leads  to  an  integral  equation  of  the  form 


/I 

fb{y)K(kx,  y,  kv)dy. 


(2.7) 


With  known  functions  Ti(kx,ky)  and  K(kx,y,ky),  Eq.  (2.7)  is  solved  numerically,  and  the  solution  gives 
both,  fi(k)  and  fu{k).  Figure  2.2(a)  is  a  contour  plot  of  the  steady-state  distribution  for  different  values 


Figure  2.2:  (Color  online)  (a)Contour  plot  of  the  distribution  functions  for  different  damping  parameters  (F  —  1 
kV/cm)  and  (b)different  applied  fields  (7  =  0.1). 


of  7  and  F  =  1  kV/cm.  For  negligible  values  of  7  (top  panel),  the  distribution  is  elongated  toward  the 
kx  axis,  with  its  centroid  located  around  kx  =  kc  ( kx  =  1  on  the  figure).  As  the  low-energy  damping 


11 


increases(7  — >  0.5),  electrons  scatter  in  the  low-energy  region  (I)  with  a  few  of  them  reaching  region  II. 
The  electron  population  concentrates  around  the  distribution  centroid  that  recedes  towards  kx  =  0  (bottom 
panel).  Notice  the  color  code  is  different  for  each  panel.  Figure  2.2(b)  display  similar  contour  plots  for 
different  electric  fields  and  constant  7  =  0.1.  Here  as  the  electric  field  is  increased,  the  distribution  is  more 
and  more  elongated,  and  broader  along  the  kx  axis,  but  an  unexpected  phenomenon  occurs  as  the  centroid 
of  the  distribution  evolves  into  a  double  hump  clearly  seen  in  the  bottom  panel.  This  is  a  consequence  of 
the  high  repopulation  rate  for  high-energy  electron  close  to  kx  =  2 kc  ( kx  =  2  on  the  figure). 

Figure  2.3  is  a  plot  of  the  current  density  in  the  high-field  regime  for  different  values  of  7.  As  expected, 
the  current  decreases  with  damping  and  reaches  it’s  saturation  value  at  higher  fields. 


Steady  State  Current  Density 


F(kV/cm) 

Figure  2.3:  (Color  online)  Current  density  in  the  high  field  regime  for  different  damping  parameters.  Here  the 
carrier  concentration  in  nc  =1.8  x  1011/cm2. 


2.4  TRANSIENT  REGIME 

The  substitution  u  =  t  —  Hkx/eF  is  used  to  solve  for  the  distribution  function  in  the  transient  regime, 

e-F  dgn{kx,ky\u)  _  gn(kx,ky;u ) 
k  dkx  rop(k) 
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where  gn(kx ,  ky\ u )  =  fu(kx,  ky,  u  +  Hkx/eF).  The  solution  of  Eq.  (2.8)  is  then 


9ii(kx,ky;u )  =  gn(k°,ky;u)exp\  -  — 


dz 


eF  Jko  TOP(z,ky)  \ 


(2.9) 


Going  back  to  the  time  variable,  the  general  solution  in  region  II  is  then 


fn(kx,kyit)  =  fn(k°,ky,t-  vQ(kx  -  fc°))exp  l 


dz 


eF  Jk o  TOP(z,ky) 


(2.10) 


where 

1  _  eF 
va  h 

is  the  “speed”  with  which  the  distribution  moves  towards  the  critical  circle.  Equation  (2.10)  describes  a 
distribution  that  starts  at  the  critical  circle  and  moves  in  the  direction  of  the  electric  field  F  with  “speed” 
l/i/Q,  while  decreasing  exponentially  as  carriers  emit  OP’s. 

At  (t  =  0),  the  distribution  in  region  I  is  the  Fermi-Dirac  distribution, 


f  i {kx  i  ky ,  t  —  0)  —  f o {kx ,  ky ) . 


(2.11) 


Right  after  the  electric  field  is  turned  on,  during  the  first  trip  towards  the  critical  circle(0  <  t  <  vakc),  the 
inside  distribution  drifts  towards  the  critical  circle  while  undergoing  low-energy  scattering.  During  this  trip, 
one  assumes  that  all  carriers  are  in  region  I,  corresponding  to  ///  ~  0.  Using  the  coordinate  transformation 
described  earlier  and  the  initial  condition,  the  distribution  in  region  I  is  given  by 


f^ikx.ky.t)  =  f0{kx  -  t/v0,  ky)e  ‘  +  —  e 

T 


/  kx  —t/vQ 


dze  °  f0(z,ky). 


(2.12) 


One  then  uses  the  boundary  condition  fi(kx,  ky,t )  =  fn(kx,  ky,t) 54  and  Eq.  (2.10)  to  obtain  the  distribution 
in  region  II: 


fn  (kx,ky,t)  =  /j1} (, k°,ky,t  -  v0 (kx  -  ) )  exp  <j  - 


dz 


kO  TOP(Z,ky)  f 


(2.13) 


The  preceding  equation  describes  carriers  in  region  II  that  interact  with  OP  and  scatter  into  region  I.  For 
t  >  vakc ,  the  solution  is  substituted  into  the  integral  of  the  right  hand  side  of  Eq.  (2.1a)  to  start 


13 


the  same  procedure,  with  the  same  variable  substitution  u  =  t  —  v0kx,  for  later  times  to  get 


fj2\kx,ky,t)  =  —e 

T 


vne 


l-k° 

pkx 


f0(z,ky)e  °  dz 


l-k° 


^2,S{k',z,ky)f[1j(k',t-  va(kx  -  z))e  °  dz 


(2.14) 


for  the  second  distribution  in  region  I.  This  distribution  is  accelerated  towards  the  critical  circle  and  the 
process  repeats  itself. 

As  the  electron  population  moves  back  and  forth  between  regions  I  and  II,  it  undergoes  each  time 
more  dephasing  and  broadening  due  to  the  finite  duration  of  the  OP  emission  process.  The  distribution 
function  in  each  region  at  any  time  t  corresponding  to  the  nth  trip  toward  the  critical  circle,  that  is,  for 
n  =  integer[eFt/hkc ]  +  1,  can  be  written  as  a  superposition  of  distributions  fjl\kx,kyit)  and  fjlj(kx,ky,t) 
of  individual  ith  “trip” ,  that  is, 

n 

fi(kx,ky,t)  =^2fjl\kx,ky,t)  (2.15a) 

2=1 

and 

n 

fii{kx,kv,t)  =  ^2fnikx,kyit).  (2.15b) 

2=1 

The  whole  distribution  is  then  normalized  to  the  total  carrier  concentration  nc. 


2.5  RESULTS 

The  dimensionless  time  parameter  t  will  be  used  to  describe  the  results: 


where 


t*  =  vakc 


is  the  approximate  time  taken  for  the  distribution  to  complete  one  “trip” .  The  results  are  presented  in  terms 
of  the  normalized  coordinates,  kx/kc  — >  kx,  ky/kc  — >  ky,  and  kc  — >  1. 

Figure  2.4(a)  shows  the  time  evolution  of  the  distribution  function  (DF)  for  F  —  1  kV/cm  (t*  «  2  ps), 
in  the  absence  of  low  energy  scattering  (7  =  0).  At  t  =  0(t  =  0),  the  initial  distribution  is  the  tail  of  the 
Fermi  distribution  in  the  conduction  band  [Eq.  (2.2)],  which  is  Maxwellian-likc,  and  drifts  along  the  kx 
direction  (opposite  to  the  F  field)  to  reach  the  OP  energy  at  t  =  1  (second  panel).  At  that  time,  a  “hump” 
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Figure  2.4:  (Color  online)  a)  3D  plots  of  the  distribution  function  at  different  times  (t  =  0,  1,  1.5,  2,  and  2.5)  for 
F  =  1  kV/cm  and  vanishing  low-  energy  scattering  (7  =  0).  b)  Corresponding  depopulation  and  repopulation  rates 
for  the  same  times  as  in  (a). 

appears  around  k  =  0  (not  yet  visible  on  the  DF  graph)  as  the  front  electrons  reach  region  II,  where  they 
emit  OP’s  and  scatter  back  inside  region  I,  as  seen  in  the  second  panel  of  Fig.  2.4(b).  As  time  progresses 
t  =  1.5,  the  “hump”  develops  into  a  dimple,  while  the  remaining  high-energy  electrons  from  the  first  trip 
continue  their  drift  in  region  II,  where  they  experience  strong  OP  depopulation.  As  seen  in  the  third  panel 
of  Fig.  2.4(b),  the  corresponding  repopulation  rate  at  low  energy  exhibits  a  craterlike  shape,  which  is  due 
to  the  randomizing  nature  of  the  deformation  potential  OP  scattering.  Indeed,  as  schematically  shown  in 
Fig.  2.1(b),  after  OP  emission  by  high-energy  electrons,  all  | k  —  kc\  values  become  equiprobable,  which 
forms  a  drifting  circle  of  high  repopulation  rate,  with  increasing  k  radius  as  first  trip  electrons  penetrate 
deeper  in  region  II.  This  craterlike  feature  of  the  repopulation  rate  is  the  primary  cause  of  the  dimple  in 
the  low-energy  distribution  function,  which  are  areas  in  k  space  where  electrons  have  low  probability  to 
scatter.  As  time  progresses,  the  DF  “dimple”  evolves  into  a  crater-like  shape,  [Fig.  2.4(a),  fourth  and  fifth 
panel] ,  and  the  successive  depopulation- repopulation  OP  processes  overlapping  at  low  energy  with  different 
amplitudes  form  also  smooth  terraces  in  the  low-energy  tail  of  the  distribution  as  it  approaches  steady  state 
[Fig.  2.4(a)  ,  fourth  panel].  These  morphological  effects  in  electron  distribution  are  unique  to  graphene 
as  a  result  its  linear  band  structure  and  the  interplay  of  the  quasiballistic  acceleration  and  relaxation  by 
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Figure  2.5:  (Color  online)  3D  plots  of  the  distribution  functions  at  t  =  2.5  for  7  =  0,  0.05,  0.1,  and  0.5.  The  applied 
field  is  F  =  1  kV /cm. 
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high-energy  monochromatic  OP’s. 


Figure  2.6:  (Color  online)  3D  plots  of  the  distribution  functions  at  t  =  2.5  for  7  =  0,  0.05,  0.1,  and  0.5.  The  applied 
field  is  F  =  1  kV /cm. 

Figure  2.5  shows  snapshots  of  the  distribution  function  at  t  =  2.5  for  varying  low  energy  relaxation¬ 
time,  expressed  in  terms  of  7  =  tq/t.  From  the  figure,  one  can  see  that  the  crater  in  the  DF  that  occurs 
for  7  =  0  progressively  disappears  as  7  increases.  Indeed  low-  energy  scattering  in  region  I  redistributes 
charge-carrier  momenta  around  k  =  0,  especially  in  the  crater  center.  For  these  reasons,  the  amplitude  of 
the  distribution  in  region  I  increases  around  k  =  0.  Also,  the  distribution  amplitude  decreases  in  region  II 
( k  >  kc )  as  carriers  spend  more  time  in  region  I  {k  <  kc).  Similarly,  the  distribution  recovers  a  streaming 
profile  along  the  kx  direction  as  conventional  semiconductors.41,54  Nevertheless,  even  for  strong  damping 
(7  =  0.5),  the  distribution  is  characterized  by  a  jagged  profile,  which  contains  the  front  and  back  ridges 
of  the  crater  remnant  still  caused  by  the  cumulative  effects  of  the  OP  scattering  for  backward  and  forward 
carrier  relaxation. 
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Figure  2.7:  (Color  online)  a)  Current  density  as  a  function  of  time  for  different  values  of  the  low-energy  scattering 
parameters  7.  The  applied  held  is  F  =  1  kV/cm.  b)  Cross  sections  of  the  distribution  at  t  =  3  and  kx  =  0.5 kc  for 
the  corresponding  values  of  7.  Dash  lines:  initial  distribution  function  at  kx  =  0. 


The  current  density  on  the  plane  is  given  by, 


Jx(t)  =  -4evF  ^2  /(&>  t)  cos(^)>  (2.16) 

k 

where  <f>  is  the  angle  between  k  and  the  kx  axis.  Figure  2.6(a)  shows  the  current  density  as  a  function  of 
time  (t)  for  F  =  1  kV/cm  and  different  low-energy  scattering  rates.  For  weak  scattering,  the  current  density 
overshoots  its  steady-state  value  through  damped  oscillations,  as  a  result  of  the  back-and-forth  motion  of  the 
distribution  function  (Fig.  2.1).  The  weaker  the  scattering,  the  higher  the  current  overshoot.  For  strong  low- 
energy  scattering,  the  current  converges  monotonically  toward  its  steady-state  value  without  oscillations.41 
Notice  that  the  stronger  the  scattering,  the  lower  the  steady-state  current  value.  Figure  2.6(b)  shows  the 
corresponding  DF  cross-sections  at  kx  =  0.5  and  t  =  3  relative  to  the  initial  distribution  at  kx  =  0.  As 
low-energy  scattering  increases,  the  DF  becomes  narrower,  taking  a  “streaming”  profile  in  the  electric-field 
direction.  This  is  due  to  the  fact  that  as  the  low-energy  scattering  increases,  fewer  electrons  reach  the  high 
energy  (  E  >  Hiuop)  region  II,  reducing  the  number  of  electrons  scattered  back  to  the  low-energy  region  I, 
which  narrows  the  distribution. 

Figure  2.7  shows  the  current  density  as  a  function  of  time  for  three  different  field  values,  with  7  =  0.05 
in  each  case.  One  observes  that  the  oscillation  period  rgr  =  Hloop / eFvF  scales  with  the  inverse  of  the 
electric  field  F.  Also  quite  expectedly,  the  overshoot  value  increases  with  electric  fields,  but  the  damping 
is  also  enhanced  with  electric  fields,  which  is  due  to  the  fact  that  the  electron  distribution  penetrates  the 
high-energy  (E  >  Hujop)  region(II)  farther  than  Hjujop,  stretching  the  carrier  relaxation  by  OP  emission, 
which  in  turn  broadens  the  DF  in  region  I  along  the  field,  thereby  reaching  steady  state  quicker.  Figure 
2.7(b)  shows  the  DF  cross  sections  at  kx  =  0.5  and  t  =  3  for  the  different  fields,  relative  to  the  initial 
distribution  at  kx  =  0.  As  the  field  is  increased,  the  distribution  function  broadens,  because  more  electrons 
reach  the  high-energy  (E  >  Hujop  )  region  II,  and  scatter  back  into  region  I  broadening  the  distribution  in 
the  process. 

Note  that,  for  Hujop  =  0.2  eV  and  F  =  1  kV/cm,  the  period  of  oscillations  rgr  =  Hujop / eFvF  =  2  ps 
corresponds  to  oscillation  frequencies  just  below  the  terahertz  range.  Also,  for  the  low  damping  case  7  =  0 
in  Fig.  2.6(a),  the  maximum  current  density  corresponds  to  carrier  velocity  v  ~  0.95vF. 

The  average  energy  density  of  carriers  in  the  conduction  band  is  given  by 

E(t)  =  4 hvF  kf(k,  t),  (2-17) 

k 
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Figure  2.8:  (Color  online)  a)  Current  density  as  a  function  of  time  for  different  fields  and  7  =  0.1.  b)  Cross  section 
of  the  distribution  function  at  t  =  3  and  kx  =  0.5 kc  for  the  corresponding  values  of  F.  Dash  lines:  initial  distribution 
function  at  kx  =  0. 
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Figure  2.9:  (Color  online)  Normalized  value  of  the  carrier  energy  as  a  function  of  time  for  different  values  of  the 
low-energy  scattering  parameters  7.  Solid  (7  =  0),  circles  (7  =  0.05),  crosses  (7  =  0.1),  and  diamonds  (7  =  0.5).  The 
applied  field  is  F  =  1  kV/cm.  b)  Same  as  (a)  but  for  different  applied  fields  and  7  =  0.1:  Solid  ( F  =  1.5  kV/cm), 
circles  (F  =  1  kV/cm),  crosses  (F  =  0.5  kV/cm). 
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where  f(k,t)  is  the  overall  normalized  distribution  function,  and  the  summation  is  taken  over  all  k  in 
regions  I  and  II.  At  t  =  0,  the  average  energy  per  charge  carrier  is  found  to  be  E{ 0)  ~  2.2/csT.  This  value 
is  roughly  twice  larger  than  ksT  expected  for  two-dimensional  systems,  and  is  the  direct  consequence  of 
the  linear  energy  dispersion  in  graphene,  by  contrast  to  the  parabolic  dispersion  in  normal  2D  systems. 
Figure  2.8(a)  displays  the  ratio  e  =  E(t)/E( 0)  for  F  =  1  kV/cm  and  different  low-energy  scattering  rate 
(7).  As  expected,  the  energy  converges  to  higher  values  as  low-energy  scattering  is  reduced.  In  addition, 
the  convergence  toward  steady  state  occurs  through  damped  oscillations,  even  for  7  =  0  as  a  consequence 
of  the  back-and-forth  motion  of  the  DF  between  the  OP  energy  and  the  carrier  zero-point  energy.  Quite 
clearly  the  oscillation  period  is  given  by  t*  for  all  7.  It  is  also  seen  that  the  oscillations  persist  even  with 
significant  low-energy  scattering.  Figure  2.8(b)  displays  the  normalized  energy  e  as  a  function  of  time  for 
different  fields  and  7  =  0.1.  As  expected,  carrier  energies  reach  higher  values  as  the  field  is  increased,  and 
the  oscillation  period  decreases  (larger  t-period). 

2.6  DISCUSSIONS 

A  transient  analysis  of  the  onset  of  current  oscillations  at  the  electric  field  turn-on  caused  by  the  back- 
and-forth  motion  of  carrier  distribution  function  between  the  zero-point  energy  and  OP’s  in  the  presence 
of  varying  damping  mechanisms  has  been  provided.  In  this  context,  the  anomalous  shape  of  the  carrier 
distribution  is  due  to  an  interplay  between  ballistic  acceleration  and  deformation  potential  OP  emission  in  the 
transient  regime.  If  OP-limited  current  oscillations  have  been  predicted  in  GaAs,41  and  indirectly  observed 
in  slightly  n-doped  InSb,42  their  manifestation  in  graphene  is  different  in  several  respects:  First,  owing  to 
the  linear  carrier  energy-momentum  dispersion,  the  oscillation  periodicity  is  given  by  Tgr  =  hujQp/eFvF 
in  graphene,  while  in  a  GaAs  parabolic  conduction  band  it  is  expressed  as  TQaAs  =  y !  2m*hjjQpAs  jeF , 
which  nevertheless  yields  similar  values,  since  the  small  effective  mass,  and  the  OP  phonon  frequency  in 
GaAs  compensate  for  the  large  vp-  41,44  Second,  in  III-V  semiconductors,  the  OP  polar  nature  focuses  the 
low-energy  repopulation  along  the  kx  axis,  which  provides  a  “streaming”  profile  to  the  carrier  distribution 
instead  of  a  crater-like  shape  in  this  case.  Finally,  in  compound  semiconductors,  the  oscillation  onset  is 
restricted  by  two  conflicting  conditions:  On  the  one  hand,  the  low  value  of  the  OP  energy  ( tkoop  ■  0.04  eV) 
is  comparable  to  the  thermal  broadening  of  the  carrier  distribution  at  room  temperature  so  that  the  back- 
and-forth  motion  of  the  distribution  between  the  optic  phonon  and  the  zero-point  energy  is  immediately 
damped. 41,43,44  On  the  other  hand,  at  low  temperature,  ionized  impurity  scattering  becomes  dominant  and 
produces  strong  damping  which  can  only  be  reduced  by  lowering  the  dopant  density,  thereby  lowering  the 
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carrier  density,  and  weakening  the  oscillation  amplitude.  In  this  respect,  the  high  conductance  of  graphene, 
and  the  high  optic  phonon  energy  provide  the  conditions  for  room-temperature  observation. 

Low-energy  scattering  should  however  still  be  minimized.  Usually,  for  experimental  studies  and  device 
applications,  graphene  layers  rest  on  a  dielectric  substrate,  or  are  confined  between  two  dielectric  slabs.55,56 
In  this  case  the  presence  of  high  K  dielectrics  sandwiching  the  strictly  2D  graphene  layer  may  be  used  to 
screen  charged  impurity  that  may  reduce  low-energy  elastic  scattering. 57  However,  dielectrics  also  contain 
interface  and  remote  static  charges  that  may  offset  dielectric  screening.58  Moreover,  the  interaction  between 
the  2D  carriers  in  graphene  and  low-energy  Remote  Interface  Phonon  (RIP)  arising  from  the  proximity  of 
the  substrate59,60  introduces  new  scattering  sources.48,61  Therefore  in  general,  suspended  graphene-layers 
avoiding  the  RIP  influence  may  be  preferable.49,62 

In  this  context  to  be  observable  at  room  temperature,  the  velocity  oscillations  should  also  take  place 
within  a  parameter  window.  On  the  one  hand  the  process  requires  eV  >  hioop  where  V  is  the  external  bias, 
to  reach  the  OP  energy,  and  on  the  other  hand,  the  field  should  be  large  enough  for  carriers  to  escape  low- 
energy  scattering,  but  huop/eFvFTOP  ^  1  where  1  /top  is  the  OP  scattering  rate,  so  that  charge  carriers 
do  not  penetrate  the  high-energy  region  E  >  Tvjjop ,  scatter  immediately  after  they  reach  the  OP  energy, 
which  maintains  the  coherence  of  the  distribution  function.  These  requirements  impose  a  lower  and  upper 
bound  on  the  electric  field,  i.e.  0.5  kV/cm  <  F  <C5  kV/cm  (in  graphene),  and  a  lower  bound  on  the  sample 
length  L  >  vftop  (>  1/rm  for  top  <  1  ps),  but  L  should  be  smaller  than  a  few  values  of  A  =  Hcoop/eF,  so 
as  to  prevent  oscillation  damping. 

One  important  issue  for  the  validity  of  this  analysis  is  the  effect  of  leakage  current  due  to  electrons  in 
the  valence  band  crossing  to  the  conduction  of  carriers  through  the  Dirac  point.  Indeed  it  has  been  shown 
that  there  is  still  a  minimum  conductance  ( G  ~  4 e2/h)  between  the  two  bands  despite  the  singular  nature  of 
the  Dirac  point.63  This  value  was  later  measured  to  be  G  ~  e2/h.24  However,  this  effect  becomes  important 
only  when  the  graphene  layer  width  w  >  23 fim,  for  which  the  band  gap  Eg  <  0.18  meV  are  vanishing. 
Moreover,  thermal  effects  such  as  acoustic  phonon  absorption  by  carriers  in  the  valence  band,  are  forbidden 
by  conservation  of  both,  energy  and  momentum.  As  for  OP  absorption  it  has  been  shown  earlier  (see  Sec. 
II)  their  occupation  number  is  also  negligible  over  the  time  scale  considered  in  this  analysis. 
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Chapter  3 


SOFT  PARAMETRIC  RESONANCE 
FOR  HOT  CARRIERS  IN 
GRAPHENE* 

3.1  INTRODUCTION 

In  high  electric  fields  F,  it  is  well  known  that  carriers  in  graphene  can  accelerate  ballistically  before  being 
scattered  by  high-energy  optical  phonons  (OP’s)  ( Hojop  ~  0.2  eV)  causing  carrier  velocity  saturation.40,39 
This  produces  a  back-and-forth  motion  of  carriers  in  k  space  between  monochromatic  OP  energy  and  the 
Dirac  point  with  a  time  period  Tp  =  Hcoop/eFvp , 41,64  which  results  in  a  carrier  velocity  overshoot65  and  even 
damped  oscillations  during  the  transient  to  steady  state,  when  the  field  is  suddenly  turned  on.41,64  While 
these  oscillations  were  predicted  to  occur  in  GaAs  at  low  temperature  so  that  hujop  A:sT[/ia;op(GaAs)  = 
36  meV]  they  were  limited  to  low  fields  (F  ~  50-100  V/cm,  tf  ~  30  ps)41  to  prevent  intervalley  scattering, 
whereas  strong  Coulomb  scattering  arising  from  the  charged  dopants  would  offset  the  effect.  In  graphene, 
the  absence  of  an  energy  gap  guarantees  carriers  without  requiring  dopants,  while  the  large  OP  energy  and 
weak  acoustical-phonon  (AP)  scattering66  allows  its  manifestation  at  room  temperature  with  ramifications 
in  THz  technology,  since  the  oscillation  frequency  wp  =  27r /rp  oc  F  ~  1  THz  (F  =  2  kV/cm)  is  tunable  with 
the  electric  field.  If  a  periodic  (ac)  field  is  superposed  onto  the  dc  field,  the  frequency  of  the  back-and-forth 
carrier  motion  is  modulated  by  the  ac  frequency,  as  a  parametric  oscillator.  As  a  result,  the  amplitude 
of  the  carrier  velocity  or  current  oscillations  is  expected  to  be  resonantly  enhanced  when  the  ac  frequency 
w  matches  a  particular  value  rj  of  the  natural  frequency  wp,  i.e. ,  uj  =  r/ujp.61  However,  there  are  distinct 
differences  between  the  usual  parametric  resonance  (PR)  and  this  type  of  hot-carrier  resonance:  First,  the 
natural  oscillations  are  strongly  damped  as  a  result  of  the  probabilistic  nature  of  the  carrier-OP  interaction 
that  relaxes  carrier  energy  at  different  times  and  momenta  once  they  reach,  and  even  overshoot  the  OP 
energy.  Second,  the  system  is  strongly  dissipative  as  the  OP  relaxation  is  responsible  for  bringing  back  the 
carriers  to  the  low-energy  Dirac  point.  Consequently,  the  resonance  is  anticipated  to  be  “soft”  i.e.  with  a 
broad  peak  in  the  oscillation  amplitude  vs  w,  and  to  manifest  for  different  values  than  normal  PR. 

Because  of  these  distinctive  features,  it  is  shown  in  this  chapter  that  the  anomalous  nature  of  this  type 

*This  chapter  closely  follows  Samwel  Sekwao,  and  Jean-Pierre  Leburton,  Phys.  Rev.  B  87,  155424  (2013). 


24 


Figure  3.1:  (Color  online)  Quasiballistic  carrier  acceleration  followed  by  OP  scattering  in  the  2D  k  space  of  graphene. 
The  circle  corresponds  to  the  electronic  energy  ftwop- 

of  resonance  that  manifests  for  77  ~  1/2  instead  of  77  ~  2  in  normal  PR.67  We  also  find  that  the  dephasing 
between  the  current  and  ac  field  exhibits  a  minimum  as  a  function  of  the  ac  field  frequency  for  weak 
damping  by  AP  or  other  low-energy  scattering,  and  softens  to  become  monotonic  at  high  damping  for  all  ac 
field  strengths. 

3.2  OPTIC-PHONON  SCATTERING  AND  HOT 
CARRIER-MODEL 

Consider  a  system  of  electrons  in  the  conduction  band  of  graphene  under  the  influence  of  a  spatially  ho¬ 
mogeneous  and  time-dependent  electric  field  F(t).  The  electric  field  takes  the  form  F(t)  =  Fa  +  F\  cos (ojt) 
and  is  applied  along  the  positive  kx  direction,  where  Fa  is  a  permanent  constant  field,  F\  is  a  constant  such 
that  0  <  F\/F0  <  1,  and  w  is  the  frequency  of  the  applied  field.  The  momentum  space  is  divided  into 
the  low  (I)  and  high  (II)  energy  regions  bounded  by  the  critical  momentum  kc  that  corresponds  to  electron 
energy  hvpkc  =  Hloop  (Fig.  3.1).  In  the  low-energy  region  I,  charge  carriers  are  ballistically  accelerated 
towards  the  critical  circle  kc  while  interacting  with  low-energy  scattering  agents  (e.g.,  AP’s  or  impurities). 
In  the  high-energy  region  II,  the  carriers  lose  all  their  energy  by  OP  emission  and  are  scattered  back  into  the 
low-energy  region.  For  this  process  to  occur,  the  electric-field  maximum  is  low  enough  such  that  electrons 
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move  back  and  forth  between  regions  I  and  II  only,  with  negligible  probability  to  reach  E  >  2 Huop- 

Because  of  the  probabilistic  nature  of  carrier  transport,  a  large-signal  Boltzmann  transport  equation 
(BTE)  that  accounts  for  low-energy  scattering  (damping),  e.g.,  by  impurities  and  AP  will  be  solved.47,64 
The  BTE  in  the  two  regions  can  be  written  as64 


dfi(k,t)  ,  eF(t)  <9/j(M)  fi{k,t)  -  fa{k)  ^  r,  r  r, 

- H7 —  +  — E - HI - = - 1-  > '  Sop(k  ,  k)fn{k  ,  f) 

at  h  okx  tle  ^ 


(3.1a) 


dfn(k,t )  eF{t)  dfn{k,t) 

dt  H  dkx 


—fn(k,  t )  ^  S0p(£,  *0) 
k' 


(3.1b) 


where  //(fc,t)  and  ///(fc,t)  are  the  time-dependent  momentum  distribution  functions  (DF)  in  the  low-  and 
high-energy  regions,  respectively,  and  e  is  the  electronic  charge.  Equation  (3.1a)  describes  electron  transport 
at  low  energy,  where  the  left-hand  side  (LHS)  accounts  for  the  transient  drift,  while  the  first  term  on  the 
right-hand  side  (RHS)  accounts  for  low  energy  scattering,  i.e,  AP  and  impurities  scattering.  The  second 
term  on  the  RHS  of  Eq.  (3.1a)  accounts  for  low  energy  carrier  repopulation  caused  OP  emission.  Equation 
(3.1b)  describes  electron  transport  at  high  energy  (E  >  hujQp),  where  the  LHS  and  the  RHS  account  for 
transient  drift,  and  electron  depopulation  due  to  OP  emission,  respectively. 64  In  Eq.  (3.1a),  the  function 
fa(k)  is  the  Fermi-Dirac  DF  (&f  >  0),  Sop(k,k')  is  the  OP  transition  rate  between  the  states  k  and  k’, 
and  tEe  is  the  relaxation  time.64  In  this  analysis,  the  temperature  of  the  graphene  sample  is  assumed  to 
be  T  =  300  K,  so  ng  <  1  and  phonon  absorption  can  be  neglected.  However,  it  is  observed  that  the  model 
is  also  valid  at  lower  T,  as  the  DF  profile  larger  than  the  thermal  broadening  is  essentially  determined  by 
high-field  carrier  dynamics,  as  long  as  the  Coulomb  scattering  (dopant  concentration)  can  be  kept  weak,  as 
shown  in  Ref.  68  and  later  on  in  this  analysis.  Note  that  at  this  temperature,  and  if  we  choose  Ep  =  kET 
above  the  Dirac  point  in  /0(fc),  the  carrier  concentration  is  nc  ~  1.8  x  1011  cm”2,  which  is  low  enough  to 
neglect  intercarrier  scattering  on  the  DF.  Moreover,  the  hole  concentration  is  even  smaller  to  significantly 
affect  the  carrier  dynamics  in  the  conduction  band  so  that  interband  transition  can  be  neglected. 69 


3.2.1  Self-Consistent  Solution  of  Boltzmann  Transport  Equation 

The  procedure  is  to  solve  Eq.  (3.1b)  for  ///(fc,t)  and  substitute  the  solution  in  Eq.  (3.1a)  to  solve  for 
fi(k,  t). 44,64  The  DFs  in  the  two  regions  are  then  matched  on  the  boundary  k  =  kc. 
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By  using  the  substitution  k  =  kx  +  j3(t),  where, 


P(t)  =  -t  j*  F(s)ds,  (3.2) 

the  LHS  of  Eq.  (3.1)  transforms  into  eF(f3~1(K—kx))dgui/Hdkx ,  where  giji(kx ,  ky\  n)  =  fiji(kx ,  ky,  /3_1(k— 
kx)),  and  /3_1  is  the  inverse  function  of  (3  so  that  /3-1/3(f)  =  t.  Consequently,  the  general  solution  of  Eq. 
(3.1b)  takes  the  form 

fii(kx,ky,t)  =  fb(ky ,/3_1(^(t)  +  kx  -  k°))M(kx,ky,t ),  (3.3) 


where  fb(ky:t)  =  fu(kx,ky,t )  is  the  time  dependent  DF  evaluated  at  the  boundary  k  =  kc,  and  kx  = 
yj ( kc )2  —  ( kv )2.  The  M(kx,ky,t)  factor  given  by 


M(kx,ky,t)  =  exp 


dp  tqp  1(p,ky)  'i 

F\p-'{p{t)+kx-p)]I 


is  the  decay  function  caused  by  OP  emission  of  hot  carriers,  and  1  /rop(fc)  =  Sop{k ,  fc').  Equation  (3.3) 

is  then  substituted  into  Eq.  (3.1a)  to  solve  for  /j(fc, t).  The  matching  conditions  fb(ky,t)  =  fi(k®,ky,t)  = 
fn{kx ,  kyi  t )  of  the  two  solutions  //  and  ///  at  the  boundary  leads  to  an  integral  equation  of  the  form 


fb(ky,t)  =  fl (ky,  t)  + 


(2ny 


l-k° 


dp  /  dk'  Sop(k',k) 


fu(k',t')ex  p{- - -}/F(t'),  (3.4) 

kx  —p  TLE 


where  the  function  fl{ky ,  t)  is  the  solution  fi{kx:  ky ,  t)  in  the  absence  of  OP  scattering,  and  t'  is  a  retarded 
time  such  that  /3(£')  =  /3(£)  +  kx  —  p  (see  Supplementary  Material  in  Ref.  70).  The  second  term  on  the 
RHS  of  Equation  (3.4)  accounts  for  the  contribution  of  OP  emission  to  the  DF  at  the  boundary,  and  the 
summation  is  taken  over  states  k'  in  the  high-energy  region.  Eq.  (3.4)  is  solved  by  iteration,  and  the  solution 
for  fb(ky,t)  is  expressed  as  a  series, 


fb{ky,t)=  fl  ( ky,t)  +  fl  (ky  ,t)  +  fl  (ky  ,£)  +  ... 


(3.5) 


which  converges  since  the  function  M(kx,  kyi  t)  is  a  decreasing  exponential  and  f{f  oc  (l/27r)"-1.  The  solution 
for  fb(ky,t)  used  throughout  this  analysis  is  obtained  by  neglecting  terms  of  0((1/27t)3)  and  higher  in  the 
series  (3.5).  Once  fb(kv,t)  is  known,  the  DFs  in  both  regions  are  readily  obtained. 
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The  two-dimensional  (2D)  current  density  on  the  plane  is  given  by 


Jx(t)  =  —4  evF  E  /(M)cos(0),  (3.6) 

k 

where  (j>  is  the  angle  between  k  and  the  kx  axis,  and  f(k,t)  is  the  DF  in  the  two  regions. 

3.3  RESULTS 

In  this  analysis,  the  value  Fa  =  1  kV/cm  is  used  and  the  applied  frequency  is  expressed  in  units  of  oj/ujf- 
A  dimensionless  damping  parameter  7  to  gauge  the  strength  of  low-energy  scattering  is  defined  as  7  = 
T~Op(k  =  1.5fcc) / r lb- 64  Also,  top  for  k  =  1.5fcc  is  chosen  as  the  intermediate  value  between  kc  and  2 kc  as 
l/r0p(fcc)  =  0. 


(a)  Fi/Fo  =  0.1;  (b)  F1/F0  =  0.8. 

Figure  3.2  shows  the  current  density  versus  time  for  different  values  of  7  and  two  field  strengths  at 
resonance,  i.e. ,  when  ui  ss  0.56wp  (see  Fig.  3.3).  It  is  seen  that  the  amplitudes  of  current  density  oscillations 
increase  as  the  applied  field  F±  increases  compared  to  F0  as  a  larger  population  of  electrons  escapes  low- 
energy  scattering  to  reach  the  OP  energy.  At  the  same  time,  electrons  also  reach  lower  velocities  during 
the  negative  cycles  of  F\.  For  this  reason,  the  current  density  swing  increases  with  Fi/F0.  One  also  notices 
distortions  in  the  current  density  oscillations  at  large  fields  [Fig.  3.2(b)]  as  the  electron  population  competes 
between  the  natural  oscillations  at  wp  and  the  oscillations  imposed  by  the  F\  field.  As  expected,  it  is 
also  seen  that  the  current  density  amplitude  decreases  with  increasing  7  as  a  result  of  increased  electron 
scattering  in  the  low-energy  region,  thereby  lowering  the  carrier  velocity. 

Figure  3.3  shows  plots  of  current  density  amplitude  versus  frequency  for  different  values  of  7  and  F\/ Fa. 
In  the  figure,  the  amplitude  is  defined  as  the  difference  between  the  maximum  and  minimum  values  of 
the  current  density.  As  seen  from  the  plots,  parametric  resonance  is  achieved  when  w/up  ~  0.56.  This 
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Figure  3.3:  (Color  online)  Current  density  amplitude  vs  frequency  for  three  different  values  of  the  damping  param¬ 
eter  7.  (a)  Ft/Fo  =  0.1;  (b)  Fi/FD  =  0.5;  (c)  Fi/F0  =  0.8. 


unexpected  result  is  due  to  the  fact  that  electrons  take  about  Tp  ~  Hkc/eF  to  reach  the  OP  energy,  and  an 
additional  Tp  to  lose  their  energy  once  they  reach  the  OP  energy,  as  they  can  still  accelerate  before  losing 
their  energy.  Consequently,  the  oscillation  period  is  about  2 Tp  and  u}resonance  ~  uip/2 ,  which  is  the  same 
as  the  current  oscillations  arising  during  the  transient  in  the  presence  of  the  dc  field  Fa  alone.64''1  This 
anomalous  value  is  due  to  the  fact  that  the  modulation  of  the  ac  field  acts  only  upon  the  first  half  of  the 
natural  period,  i.e.,  when  carriers  are  field  driven  toward  the  OP  energy,  while  the  second  half  of  the  period 
when  OP’s  are  emitted  is  stochastic  with  a  more  complicated  dependence  on  the  field  [see  Eq.  (3.3)],  which 
is  why  the  PR  frequency  is  not  exactly  half  of  the  natural  frequency  tup.  From  the  figure,  the  oscillation 
amplitudes  increase  with  7,  which  as  explained  in  Fig.  3.2  is  due  to  increased  scattering  at  low  energy, 
sending  back  electrons  close  to  the  k  =  0  region,  thereby  further  reducing  the  minimum  values  of  the  current 
density.  The  maximum  values  of  the  current  density  are  not  as  affected  because  a  substantial  population  of 
electrons  is  still  able  to  reach  high  energies,  even  at  high  7.  Also,  it  is  seen  from  the  plots  that  the  amplitude 
increases  with  F\/ F0 ,  as  expected,  since  the  difference  between  current  density  maxima  and  minima  increases 
with  F\/F0.  Even  though  PR  is  achieved,  it  is  rather  “soft”  because  of  the  strongly  dissipative  nature  of 
the  back-and-forth  motion  of  charge  carriers  in  the  constant  field  followed  by  OP  emission.  Obviously,  this 
effect  is  more  pronounced  for  the  higher  values  of  7  (low-energy  scattering)  and  F\/F0  (OP  scattering)  seen 
in  the  figure. 

Figure  3.4  shows  normalized  current  densities  and  electric  fields  versus  time  for  different  values  of  the 
parameters  uj/ujf  and  7  at  low  electric  fields  ( F\/Fa  =  0.1,  left  column),  and  high  fields  ( F\/Fa  =  0.8, 
right  column).  At  low  fields,  the  current  is  sinusoidal  as  expected  from  the  linear  response  to  the  field.  It 
is  also  observed  that  at  very  low  frequencies  ( oj/ujf  =  0.001)  and  in  the  quasiballistic  regime  (7  =  0.01, 
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Figure  3.4:  (Color  online)  Normalized  values  of  the  oscillating  part  of  the  current  density  vs  time.  Left  column, 
F\/F0  =  0.1;  right  column  F\/F0  =  0.8.  Top  row:  7  =  0.01;  middle  row:  7  =  0.1;  and  bottom  row:  7  =  0.3.  Blue 
curve:  uj/ujf  =  0.001;  red  curve:  lo/ujf  ~  0.56;  green  curve:  lo/ujf  =  2.3;  and  black  curve:  normalized  ac  field. 
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Figure  3.5:  (Color  online)  2D  color  plot  of  the  electron  distribution  function  difference  A /  (as  defined  in  the  text, 
and  normalized  to  the  carriers  density)  in  k  space  (normalized  units  of  k/kc)  at  four  different  times  for  field  ratio 
F\/Fa  =  0.1  and  damping  7  =  0.01.  Dashed  circles  correspond  to  the  boundary  k  —  kc  between  the  low-  and 
high-energy  regions.  Top  row:  uj/ujf  =  0.001.  Bottom  row:  uj/lof  =  0.56. 
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top  left  column),  the  current  density  is  180°  out  of  phase  with  the  field.  In  this  case,  as  the  field  slowly 
decreases  from  i  =  0toi  =  7r/w,  the  electronic  system  evolves  adiabatically  from  a  regime  in  high  fields  to 
that  in  low  fields,  which  depletes  the  charge  carriers  in  the  high-energy  region  ( k  >  kc )  and  increases  their 
concentration  in  the  low-energy  region  ( k  <  kc).  The  current  increases  as  the  number  of  electrons  with  high 
kx  values  (kc/ 2  <  kx  <  kc )  in  the  low-energy  region  increases  as  a  result  of  quasiballistic  transport  that 
results  in  a  streamed  DF.  This  situation  is  clearly  seen  in  Fig.  3.5  (top  row,  first  two  panels),  which  shows 
the  change  in  the  DF  with  time  A /  =  f[uit  =  (n  +  1)7t/2]  —  /[cut  =  mr/2],  From  t  =  f/u>  to  t  =  2i t/u), 
the  current  decreases  as  the  field  increases  because  the  electrons  that  penetrate  deep  into  the  high-energy 
region  ( k  >  kc )  with  kx  ky  are  scattered  by  OP  emission  equally  to  all  k'  =  k  —  lo/vf  values.  Indeed,  the 
absence  of  q  (phonon  wave-vector)  dependence  in  the  deformation  potential  OP  matrix  element  contributes 
to  randomizing  the  DF,50  specifically  populating  low-energy  k  states  away  from  the  field  direction.  This 
effect  results  in  lowering  the  current  (Fig.  3.5,  top  row,  last  two  panels).  In  Fig.  3.4  (left  column,  top), 
it  is  seen  that  the  phase  between  the  current  and  the  field  reaches  a  minimum  for  frequencies  approaching 
resonance  (oj/uip  =  0.56).  This  effect  is  better  understood  as  the  field  increases  from  t  =  f/uj  to  t  =  3tt/2oj, 
then  to  t  =  2i r/w,  when  the  ac  and  the  dc  fields  combine  to  enhance  the  back- and- forth  motion  of  carriers 
between  the  low-energy  k  states  (Fig.  3.5,  bottom  row,  third  panel)  and  the  high-energy  k  states  (fourth 
panel  and  first  panel).  For  frequencies  higher  than  resonance  (oj/ujf  0.56),  the  dephasing  between  current 
and  electric  field  starts  to  increase  again  (Fig.  3.6). 

One  also  observes  that  as  low-energy  damping  increases  (7  =  0.1  and  7  =  0.3;  Fig.  3.4,  left  column), 
the  dephasing  between  the  current  and  the  field  at  frequencies  below  resonance  decreases  (Fig.  3.6),  as 
low-energy  collisions  result  in  diffusive  transport  that  scatter  electrons  with  high  kx  states,  thereby  changing 
the  streamed  DF  into  a  wider  ( ky  states)  DF  with  lower  current  density.  At  intermediate  damping  (7  =  0.1), 
there  is  even  a  slight  maximum  before  resonance,  but  above  resonance  the  dephasing  increases  monotonically 
for  all  damping. 

In  higher  ac  fields  (Fig.  3.4,  right  column),  aside  from  the  fact  that  the  current  curves  are  distorted  by 
transport  nonlinearity  caused  by  competition  between  the  dc  and  ac  fields,  the  results  concerning  the  phase 
difference  between  the  oscillating  F\  field  and  the  current  densities  are  qualitatively  the  same.  One  notices, 
however,  that  the  distortions  do  not  affect  the  current  at  resonance,  which  remains  quasisinusoidal.  The 
effect  of  7  on  the  current  density  phase  is  also  seen  in  Fig.  3.6.  As  expected,  current  density  lags  behind 
the  electric  field  before  resonance,  but  the  dephasing  also  drops  around  resonance  for  low  damping. 


31 


3.5 


3 
2.5 
2 

®  1.5 

1 
0.5 
0 

-0.5 

0  0.5  1  1.5  2  2.5 

co/coF 

Figure  3.6:  (Color  online)  Phase  of  the  current  density  [compared  to  F(t)]  for  different  values  of  the  damping  7 
and  field  ratio  F\/F0  =  0.1. 

3.4  CONCLUSIONS 

Although  this  analysis  is  performed  for  time-dependent  ac  fields  in  the  condition  of  spatial  uniformity,  it  is 
also  valid  in  the  inverse  condition  of  long-range  periodically  (oscillatory)  modulated  potential  V ( x )  =  V (x+d) 
in  the  steady  state.  This  can  be  seen  from  Eqs.  (3.1),  where  the  time-dependent  differential  d/dt  operator 
of  the  BTE  LHS  is  replaced  by  the  spatially  varying  vf  cos  (f>d/dx  operator,  for  which  (j>  ~  0  in  streamed 
DFs.  Therefore,  by  making  the  substitution  t  —1  x/vf  in  the  formalism,  the  resonance  condition  between 
the  periodic  potential  and  the  hot-carrier  dynamics  will  arise  for  F„  =  (fouop)/0.56ed  in  the  presence  of  an 
external  field  Fa,  which  could  be  used  as  field  detector. 

Let  us  notice  that  for  carrier  oscillations  to  occur  in  graphene,  the  electric  field  has  to  be  high  enough  to 
escape  low-energy  scattering  agents,  but  not  too  high  as  to  overshoot  the  OP  energy. 64  Also,  at  high  damping, 
more  scattering  occurs  in  the  low-energy  region  and  the  resulting  oscillations  are  promptly  damped.  This 
problem  persists  even  with  the  addition  of  a  high  amplitude  ac  field,  and  hence  the  wider  resonance  peak 
at  high  7  in  Fig.  3.3.  As  a  result,  clean  graphene  samples  should  be  used  with  low  values  of  F\/F0  for  this 
effect  to  be  observed,  and  be  the  basis  for  novel  device  applications  either  as  a  THz  source  or  detector. 


y  =  0.01 
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Chapter  4 


ELECTRICAL  TUN  ABILITY  OF 
SOFT  PARAMETRIC  RESONANCE 
BY  HOT  ELECTRONS  IN 
GRAPHENE* 


Recently,  it  has  been  suggested  that  the  electron  current  may  exhibit  damped  transient  oscillations  due  to 
the  sudden  quasi-ballistic  acceleration  of  electrons  up  to  the  optic  phonon  (OP)  energy  Swop,  and  their 
subsequent  relaxation  by  OP,  when  a  DC  electric  field  Fa  is  turned-on.64  This  effect  would  arise  provided 
that  low-energy  scattering  is  weak  and  the  field  is  not  too  high  i.e.  Fc  ~  0.5-2  kV/cm,  yielding  oscillation 
frequencies  of  the  order  of  uip  =  eF0vp /hioop  in  the  terahertz  range.64  It  was  also  found  that  when  an  ac 
field  is  superimposed  to  the  dc  field  a  resonance  in  the  oscillations  occurs  when  the  frequency  of  the  ac  field 
ui  =  r]iOF  with  r]  =  0.56  for  a  field  F0  =  1  kV /cm.  This  anomalous  effect  was  called  soft  parametric  resonance 
(SPR)  because  the  frequency  of  the  natural  oscillations  is  modulated  by  the  ac  field  in  a  strongly  dissipative 
electronic  system.68  Asides  from  its  fundamental  peculiarity  where  the  SPR  coefficient  r j  is  anomalously 
smaller  than  one,  unlike  in  normal  PR  where  it  is  exactly  2,6'  SPR  offers  some  technological  advantage,  as 
u}f  is  tunable  with  the  DC  field  Fa.  In  this  chapter  an  investigation  the  sensitivity  the  SPR  coefficient  77  to 
the  DC  electric  field  and  the  temperature  for  possible  applications  in  THz  technology  is  provided.  For  this 
purpose,  the  linear  response  theory  that  is  shown  to  be  less  tedious  than  the  large  signal  analysis  previous 
used  to  predict  SPR,  but  is  also  in  good  agreement  with  it  is  used.  This  model  not  only  shows  a  weak 
sensitivity  of  77  to  F0  and  temperature,  but  also  reveals  a  new,  but  smaller  intensity  resonance  w  =  lof,  i.e., 
77  =  1  at  lower  T. 

Consider  a  system  of  electrons  in  the  conduction  band  of  graphene,  under  the  influence  of  a  spatially 
homogeneous  and  time  dependent  electric  field  F(t).M  The  electric  field  is  applied  in  the  positive  kx  direction 
and  is  of  the  form  F(t)  =  F0  +  ,  where  fq  is  the  ac  amplitude  with  fq  <C  Fa.  As  in  previous  work,64,68 

define  two  regions  i.e.  a  low  (/),  and  the  high  (II)  energy  regions  in  the  momentum  space,  separated  by 
the  critical  momentum  kc  =  ujop/vf  (Refs.  64  and  44)  corresponding  to  the  electron  energy  Hojop  =  hvpkc 
above  which  electrons  emit  OP’s,  and  scatter  back  to  region  I ,  where  they  undergo  quasi-ballistic  acceleration 

*This  chapter  closely  follows  Samwel  Sekwao,  and  Jean-Pierre  Leburton,  Appl.  Phys.  Lett.  103,  143108  (2013). 
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and  low  energy  scattering  with  acoustic  phonon  and  impurities. 

The  relaxation  time  approximation  is  used  to  account  for  low  energy  scattering  (e.g.,  impurities  and 
A. P), 47,64  so  the  Boltzmann  Transport  Equation  (BTE)  in  the  two  regions  read,64 


dfi{k,t )  eF(t)  df!(k,t) 

dt  h  dkx 


fi(k,t)  -  fFp{k) 

tle 


+  ^2  Sop(k' ,  k)fn(k',t), 

k' 


(4.1a) 


dfujk,  t)  eF(t)  dfn(k,t) 

dt  H  dkx 


—fu(k,  t)  ^2  S0p(k,  k'), 

k' 


(4.1b) 


where  fi{k,t)  and  fn(k,t)  are  the  distribution  functions  in  the  low  and  high  energy  regions,  respectively. 
The  function  /pp>(/c)  is  the  Fermi-Dirac  distribution  function  ( kp  >  0),  Sop{k:k')  is  the  OP  transition 
rate  between  the  states  k  and  k' ,  and  tle  is  the  relaxation  time.64  As  in  previous  work,  OP  absorption  is 
neglected  as  Nq  -C  1  (kpr  <C  Hloop )•  Since  F0  -C  Fi,  we  look  for  solutions  f(k,t)  =  f°(k)  +  f1(k)elult ,  where 
f°(k)  is  the  steady  state  solution  of  Eq.  (4.1)  with  the  applied  field  E0,64  and  /7(fc)  f°(k).  Substituting 

this  expression  into  Eq.  (4.1),  and  linearize  the  BTE  we  get  the  following  equations  for  /1(fc), 

+  ^50p(fc',fc)//1/(fc'),  (4.2a) 

k' 


^  Sop(k,  k')  +  iuj 
k' 

where  f}(k)  and  fjj(k)  are  the  solutions  fl(k)  in  the  low  and  high  energy  regions,  respectively.  The  functions 
ff{k)  and  ffj(k)  are  the  steady  state  distribution  functions  in  the  low  and  high  energy  regions,  respectively. 
Equation  (4.2b)  is  readily  solved,  and  the  solution  f}j(k )  is  then  substituted  into  Eq.  (4.2a)  to  solve  for 
f}(k).  The  two  solutions  are  then  matched  at  the  boundary  k  =  kc  to  obtain  a  solution  everywhere  on  the 
plane.44,64,68  Then  we  compute  the  2D  current  density  on  the  plane  given  by 


eFa  df}j(k) 
H  dkx 


)  (p)  eF i  dfu(k) 

)tnW-—pr  dkx  ’ 


(4.2b) 


eK9_m + ( j_ +iu)I]{i) = 

n  akx  TLE  n  okx 


Jx{t)  =  —  4ei>F  ^  5i/(fc,  t)cos(^),  (4.3) 

k 

where  <j>  is  the  angle  between  k  and  the  kx  axis,  and  f(k,t)  is  the  distribution  functions  in  the  two  regions. 
In  this  analysis,  the  results  are  expressed  in  the  normalized  ac  frequency,  uj/ojf,  where  wp  =  and 


tf 


frlOOP 

eFvF 


(4.4) 
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is  the  time  taken  by  ballistic  carriers  to  reach  the  OP  energy  tkuop  from  the  Dirac  point.  Also,  a  dimensionless 
damping  parameter  7  given  by 

7  = - ^  -  (4.5) 

tle  Sop{k,  fc')|fc=i.5fec 

is  used  to  gauge  the  strength  of  scattering  in  the  low  energy  region.64,68  Fig.  4.1(a)  shows  plots  of  the  steady 


k 

X 

Figure  4.1:  (Color  online)  (a)  Cross-sections  of  the  steady  state  distribution  function  f°  for  different  values  of 
Fa.  7  =  0.01.  Solid  lines  ( Fa  =  0.5  kV/cm),  dashed  lines  ( Fa  =  1  kV/cm),  dotted  lines  ( Fa  =  1.5kV/cm),  and 
dashed-dotted  lines  ( F0  =  2.0  kV/cm).  (b)  Real  part  of  the  solution  / 1  for  different  values  of  uo/ljf ■  Solid  lines 
(lo/ujf  =  0.2),  dashed  lines  (uj/ujf  =  0.54),  dotted  lines  (oj/ujf  =  1.08),  and  dashed-dotted  lines  {uj/ujf  =  1.5). 
7  =  0.01.  (b)  Imaginary  part  of  the  solution  f1  for  different  values  of  uj/ujf ■  Solid  lines  (cj/cvf  =  0.2),  dashed  lines 
( uj/ujf  =  0.54),  dotted  lines  (uj/lof  =  1.08),  and  dashed-dotted  lines  (c ij/lof  =  1-5).  7  =  0.01. 


state  distribution  f°  along  the  kx  axis,  for  different  values  of  the  applied  field  F0,  and  7  =  0.01.  As  seen  from 
the  figure,  the  concentration  of  carriers  in  region  II  increases  as  F0  increases.  This  is  because  carriers  are 
more  likely  to  penetrate  deeper  into  the  high  energy  region  before  they  are  scattered  by  OP’s.  Figures  4.1(b) 
and  4.1(c)  show  plots  of  the  real  and  imaginary  parts  of  the  solution  / 1  respectively,  for  different  values  of 
the  applied  frequency  uj/ujf,  and  7  =  0.01.  One  notices  the  real  and  imaginary  parts  of  / 1  oscillates  as  a 
function  of  kx  with  decreasing  amplitudes  as  the  applied  frequency  lo/ujf  increases.68  We  believe  the  origin 
of  these  oscillations  resides  in  the  coupling  between  time  and  fcx-momentum  during  the  drift-motion  of  the 
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distribution  in  /e-space,  reflected  in  the  “pseudo” -time  invariance  transformation  k'x  =>  kx  +  f*  F(s)ds  of  the 
differential  terms  of  the  BTE  in  Ref.  68 


y  =  0.01  y  =  0.1  y  =  0-3 


Figure  4.2:  (Color  online)  Current  density  amplitude  versus  applied  frequency  for  different  values  of  7.  Comparing 
results  from  the  Large  Signal  Analysis  and  the  Small  Signal  Analysis  at  T  =  300  K.  Also  shown  SSA  results  at 
T  =  77  K.  F0  =  1  kV/cm,  and  Fx/F0  =  0.1.  (a)  7  =  0.01;  (b)  7  =  0.1;  (c)  7  =  0.3. 

Fig.  4.2  shows  a  comparison  of  current  density  amplitude  versus  applied  frequency  results  obtained  by 
full  solution  of  the  BTE,  and  referred  here  as  Large  Signal  Analysis  (LSA),68  and  the  linear  solution  obtained 
from  Eq.  (4.2),  referred  as  the  Small  Signal  Analysis  (SSA),  at  T  =  300  K.  The  figure  also  shows  SSA  results 
at  T  =  77  K.  For  both  temperatures,  the  applied  field  is  Fa  =  1  kV/cm,  and  F\/Fa  =  0.1.  From  the  figure, 
one  can  see  that  the  two  methods  (LSA  &SSA)  give  virtually  the  same  results  for  the  7-values  considered, 
as  the  two  curves  are  on  top  of  one  another,  thereby  confirming  the  validity  of  the  linear  approximation  for 
F\/F0  =  0.1.  One  notices  the  SPR  at  T  =  77  K  persists  at  u/ujf  ~  0.56, 68  but  is  a  little  narrower  than  at 
T  =  300  K,  while  a  secondary  peak  is  developing  at  co/wf  ~  1.  This  second  peak  is  due  to  the  resonance 
of  carriers  oscillating  between  the  Dirac  point  and  the  OP  energy  fkoop  for  which  the  frequency  is  given 
by  up-  At  low  temperatures,  the  carrier  distribution  function  is  narrower,  which  makes  carrier  acceleration 
toward  OP  more  coherent.  One  notices  however,  that  this  second  resonance  does  not  involve  the  time  spent 
by  electrons  to  emit  OP’s  at  high  energy,  but  only  the  acceleration  toward  OP  energy.  In  fact,  it  is  also 
present  as  a  weak  shoulder  at  T  =  300  K,  but  is  revealed  at  T  =  77  K  as  the  SPR  sharpens  due  to  the 
narrower  distribution.  Note  that  the  LSA  and  SSA  results  at  T  =  300  K  are  obtained  for  frequencies  in  the 
range  0.1  <  uj/uif  <  2,  while  SSA  results  at  T  =  77  K  are  for  the  frequencies  0  <  ui/ujf  <  2. 

Fig.  4.3  shows  the  SSA  current  density  amplitude  versus  applied  frequency  for  different  F0  values  at 
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Figure  4.3:  (Color  online)  Current  density  amplitude  vs  applied  frequency  for  different  values  of  F0.  (a)  T  =  77  K. 
(b)  T  =  300  K.  F\/F0  =  0.1,  and  7  =  0.01.  Solid  line  ( Fa  =  0.5  kV/cm),  dashed  lines  ( F0  =  1  kV/cm),  dotted  lines 
(F0  =  1.5  kV/cm),  and  dashed-dotted  lines  ( F0  =  2  kV/cm). 

T  =  77  K  and  T  =  300  K,  and  for  F\jF0  =  0.1,  and  7  =  0.01.  For  both  temperatures  the  SPR  peak  shift 
slightly  to  lower  frequency,  passing  by  a  small  maximum  at  F„  =  1  kV/cm  as  Fa  increases,  thereby  indicating 
that  the  rj  parameter  is  a  weak  function  of  the  DC  field,  as  discussed  in  Fig.  4.4.  Also  for  both  temperatures, 
the  second  resonance  peak,  while  following  the  same  trend  than  the  SPR  peak  is  more  prominent  at  low 
field. 

In  Figure  4.4,  left  axis,  it  is  shown  the  variation  of  both  resonance  frequencies  versus  applied  DC  fields 
for  both  temperatures  for  which  both  resonance  peaks  remain  practically  unchanged.  It  can  be  seen  that 
the  SPR  frequency  (77  parameter)  decreases  with  the  applied  F„  field  for  both  temperatures  .  As  shown  in  a 
previous  work,68  SPR  is  achieved  when  the  applied  frequency  to  is  close  to  the  actual  frequency  of  damped 
oscillations  of  the  carrier  distribution  during  the  transient  turn-on  of  the  DC  field.68  The  actual  period 
t*  of  these  oscillations  consists  of  the  quasi-ballistic  carrier  acceleration  in  the  low  energy  region,  tf,  and 
the  time  spent  in  the  high-energy  region  before  emitting  OP’s.  As  seen  from  Fig.  4.1(a),  as  Fa  increases, 
electrons  are  likely  to  spend  more  time  in  the  high  energy  region,  before  OP  emission  occurs,  and  the  period 
t* /rp  also  increases.64  As  a  result,  one  expects  the  resulting  SPR  frequency  to  decrease  with  Fa.  However 
this  variation  of  the  parameter  is  relatively  weak  as  seen  in  Figures  4.3  and  4.4,  which  is  important  for  the 
SPR  applications  in  terahertz  technology,  since  it  indicates  that  the  SPR  frequencies  are  tunable  with  F0  at 
both,  low  and  room  temperatures.  The  right  axis  of  Fig.  4.4  shows  the  second  resonant  frequency  u>2  versus 
applied  fields,  which  is  practically  constant  at  W2 /wf  =  1,  except  at  very  low  fields,  where  it  exceeds  ujf- 
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Figure  4.4:  (Color  online)  Resonant  frequency  versus  applied  field  for  both  T  =  77  K,  and  T  =  300  K.  Left  axis: 
First  peak  (solid).  Right  axis:  Second  peak  (dash).  7  =  0.01. 


Figure  4.5  displays  current  density  versus  applied  frequency  for  different  values  of  7,  at  Fa  =  0.5  kV/cm 
for  both  temperatures.  In  both  figures  peaks  are  sharper  at  T  =  77  K  than  at  T  =  300  K,  and  as  7 
increases,  the  first  resonance  peak  shifts  to  lower  values  of  cuspr/ojf,  as  stronger  low  energy  scattering 
delays  the  carrier  distribution  in  reaching  the  high  energy  region  to  emit  OP’s,  thereby  decreasing  the 
resonant  frequency  ujspr/uf,  as  already  explained  in  Ref.  68.  As  in  Figs.  4.2  and  4.3,  the  second  resonance 
peak  is  also  more  pronounced  at  T  =  77  K  than  at  T  =  300  K.  However,  unlike  the  SPR  peak,  the  second 
resonance  peak  shifts  to  higher  values  of  uj/ujf  as  7  increases  from  0.01  to  0.1. (  Fig.  4.5(a),).  Here  too,  strong 
scattering  at  low  energy  slows  down  carriers  to  reach  the  OP  energy  at  a  higher  rate,  thereby  increasing 
u:/ojf  for  resonance.  The  same  effect  is  seen  in  Fig.  4.5(b)  where  the  shoulder  that  identifies  the  second 
resonance  shifts  towards  higher  values  of  uj/luf  as  7  increases.  At  very  high  values  of  7(7  =  0.3),  and  both 
temperatures  the  oscillating  current  amplitudes  increase,  mainly  due  to  the  fact  that  during  the  low  cycle  of 
the  amplitude  the  negative  ac  field  reduces  the  overall  field,  which  enhance  the  effect  of  low  energy  scattering 
as  already  observed  in  Ref.  68.  In  addition,  carrier  acceleration  in  the  DC  field  Fa  struggles  to  overcome 
low  energy  scattering,  so  the  second  resonance  is  less  pronounced. 

It  has  been  shown  that  parametric  resonance  is  achieved  at  different  frequencies  for  different  values  of 
the  applied  dc  field  F0,  and  ac  fields  with  amplitude  F\/Fa  =  0.1.  The  resonance  is  thus  tunable  with  F0  at 
low  and  room  temperatures.  These  results  could  have  potential  applications  in  terahertz  technology.  The 
emergence  of  a  second  resonance  peak  at  low  temperatures  is  also  observed.  This  resonance  peak  is  due  to 
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Figure  4.5:  (Color  online)  Current  density  amplitude  vs  frequency  for  different  values  of  7.  (a)  T  =  77  K  and  (b) 
T  =  300  K.  Fi/F0  =  0.1,  and  F0  =  0.5  kV/cm.  Solid  lines  (7  =  0.01),  dashed  lines  (7  =  0.1),  and  dotted  lines 
(7  =  0.3). 


the  acceleration  of  carriers  followed  by  OP  relaxation  at  energy  Hu>op-  If  Fa  is  too  high,  the  system  may 
not  reach  resonance  because  of  thermal  broadening  and  or  carriers  escaping  the  high  energy  region  without 
OP  scattering.  For  best  results,  cleaner  graphene  samples  should  be  used  with  the  appropriate  dc  field. 
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Chapter  5 


TERAHERTZ  HARMONIC 
GENERATION  IN  GRAPHENE* 


An  interesting  effect  could  arise  if  carriers  in  graphene  are  placed  in  periodic  long  range  and  time  varying 
scattering  to  achieve  transport  resonance  as  well  as  possible  frequency  mixing.  This  kind  of  situation  could 
be  realized  in  free  standing  graphene  sheets  lying  over  periodically  spaced  narrow  electric  gates  that  would  be 
regulated  by  an  a-c  field  of  appropriate  frequency  to  modulate  coulomb  scattering  of  remote  oxide  impurities 
when  carriers  pass  in  front  of  the  gates  (Fig.  5.1).  In  the  device,  a  DC  field  Fa  is  set  up  between  the  source  S 
and  the  drain  D ,  and  the  a-c  field  Fi  is  applied  between  successive  gates,  so  that  charge  carriers  experience 
periodic  electric  fields  and  scattering  times  varying  in  time  and  distance  along  the  channel. 

In  this  chapter,  it  is  shown  that  electronic  current  in  graphene  under  the  influence  of  time  and  space 
dependent  periodic  scattering  and  electric  fields  exhibits  sharp  resonances  in  the  terahertz  range,  associated 
with  the  generation  of  higher  harmonics  characterized  by  large  Q-factors  tunable  with  scattering  periodic¬ 
ity.  The  electron  dynamics  is  investigated  with  a  semi-classical  Boltzmann  formalism,  where  the  effects  of 
electron-electron  interactions  are  ignored  for  low  carrier  concentrations  (n  <C  1012/cm2).2'  Also,  the  electric 
fields  considered  in  this  study  are  not  strong  enough  so  that  interband  effects  due  to  impact  ionization72 
and  tunneling73  can  be  ignored. 

As  in  the  previous  chapters,  consider  electrons  in  the  conduction  band  of  graphene  in  a  field  of  the 
form  F(x,t)  =  F0  +  F\e with  0  <  F\/Fa  <  1,  where  Fa,  and  F\  are  the  DC  and  a-c  components 
respectively.  The  parameter  q  is  the  electric  field  wave  number,  and  u  is  the  a-c  frequency.  In  such  fields, 
electrons  accelerate  until  they  gain  enough  energy  to  interact  with  OP’s  once  E  >  Hujop,  and  suddenly  lose 
their  energy  and  momentum  to  re-accelerate  quasi-ballistically  in  the  fields,  as  this  process  repeats  itself.  In 
this  dynamical  picture,  the  momentum  space  can  be  divided  into  two  regions  /,  and  II  that  corresponds  to 
electron  energies  E  <  Ed  op  ,  and  E  >  Ejjop  respectively.  In  region  J,  electrons  accelerate  while  interacting 
with  low  energy  scattering  agents  (AP’s  and  impurities).  In  region  II,  the  electrons  interact  with  OP’s  and 
scatter  back  to  region  I.  In  low  electron  concentration  (n  <  1012/cm2),  most  of  the  electron  population  is 
stretched  along  the  fields  (streaming  case),  so  that  the  majority  of  carriers  have  their  velocity  pointing  in 

*This  chapter  closely  follows  Samwel  Sekwao,  and  Jean-Pierre  Leburton,  Appl.  Phys.  Lett.  106,  063109  (2015). 
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Figure  5.1:  Color  online)  Schematics  of  a  graphene  based  FET  with  periodic  gating.  The  DC  field  F0  is  applied 
between  the  source  and  the  drain,  and  the  AC  field  F\  is  applied  between  successive  gates.  The  wavelength  of  the 
AC  field  is  twice  the  distance  between  successive  gates. 

that  direction.64,73’74,75  If  the  fields  are  applied  along  the  x-  direction,  the  time  and  space  dependent  BTE 
in  both  regions  reads 


Ofi  df j  eF(x,  t)  df / 

dt  F  dx  h  dkx 


^  +  ^2S0p(k,,k)fII{k',x,t ) 

TyE-,  t) 

k' 


(5.1a) 


dfn  ;  dfu  eF(x,t)  dfu 
dt  VF  dx  h  dkx 


-fn  ^2  Sop{k ,  k') 
i v 


(5.1b) 


where  //,  and  fjj  are  the  distribution  functions  in  regions  I  and  II  respectively,  and  fpjj  is  the  Fermi- 
Dirac  equilibrium  distribution  function  ( kp  >  0).  Also,  in  the  streaming  case,  one  can  assume  vx  ~  vf- 
Low  energy  elastic  collisions  in  region  I  with  impurities  and  AP’s  are  accounted  for  within  a  local  and  time 
dependent  relaxation  time  approximation  64r (x,  t)  with  the  same  periodicity  as  the  gate  spacing  and  applied 
a-c  field  i.e. 

— - -  =  —  +  — ei(gx-wt)  (5.2) 

r(x,  t)  T0  n 

where  tq  is  a  constant  relaxation  time  used  to  modulate  the  strength  of  low  energy  scattering  (with  impurities 
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and  AP’s),  and  n  modulates  scattering  due  to  the  periodic  gates  (see  Fig  5.1).  The  expression  Sop(k,k') 
is  the  OP  transition  rate  between  the  states  k  and  k'  given  by52 


Sop(k,k') 


nD02(Nq  +  1)(1  +  cos(0')) 
aAuiop 


S(E ' 


E  +  fiioop) 


(5.3) 


where  Da  is  the  optical  deformation  potential,  Nq  is  the  phonon  occupation  number,  6 1  is  the  angle  between  k 
and  k',  a  is  the  area  density  of  the  graphene  sheet,  and  A  is  the  area  of  the  sheet.  The  effects  of  optic  phonon 
absorption  are  neglected  since  their  population  is  negligible,  with  Hluop  kpT  even  at  room  temperature. 
The  second  term  on  the  Right  Hand  Side  (RHS)  of  Eq.  5.1a  is  due  to  carriers  repopulating  region  I  after 
scattering  with  OP’s,  and  the  RHS  of  Eq.  5.1b  is  the  corresponding  depopulation  term. 

Due  to  the  nature  of  the  problem,  a  solution  of  the  form 


OO 

f(k,x,t)  =  fh{k)  + 


^  ^  fmn 


fmn(k)e 


i(mqx—nojt ) 


m——oo  n=— oo 


(5.4) 


is  assumed,  where  fh(k)  =  /°°  is  the  solution  to  the  homogeneous,  steady  state  problem,64  fmn(k )  are 
the  harmonic  amplitudes,  and  the  summation  excludes  all  terms  with  m  =  0  and  n  =  0,  since  they  would 
result  either  in  homogenous  (to  =  0)  or  state-state  (n  =  0)  terms  that  are  already  taken  into  account 
in  fh(k)  =  f00(k).  Substituting  Eq.  5.4  into  Eq.  5.1,  leads  to  the  following  recurrent  equations  for  the 
individual  harmonics; 


eF0  Of?171 
H  dkx 


- b  i{mqvF 

-  7~o 
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eFi  jjrn-l)(n-l) 

H  dkx  n 

'Z,sOP(k',k)f?in(k ') 

k' 


(5.5a) 


eFa  dfY}n 
h  dkx 


S'op(fc,  k ')  +  i(mqvF  —  nui) 
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eF1  df^-1){n-1] 
fi  dkx 


(5.5b) 


where  /j71",  and  fjjn  are  the  harmonic  amplitudes  in  regions  I  and  II  respectively.  On  setting  m  =  n  =  0 
into  Eq.  5.5,  the  terms  in  /°°  satisfy  the  steady  state,  homogenous  Boltzmann  equation  in  the  two  regions, 
so  that  we  get  the  following  equations  for 


efi  d/71-1  f-1'-1 

h  dkx  t\ 


(5.6a) 
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eFx  dfJT L’~1 
h  dkx 


(5.6b) 


=  0. 


Eq.  5.6b  has  the  general  solution  fjj1,  (kx,k y)  =  g(ky )  which  is  unphysical,  and  should  be  set  to 
zero,  which  by  continuity  in  region  I  also  yields  fj  =  0.  Because  of  the  recurrence  between  fmn  and 
jn  gq.  (5.5),  we  also  have  fmm  =  0  for  all  m  <  —  1. 

Let  us  now  consider  the  harmonics  fln  with  n  >  1.  Since  the  R.H.S  of  Eq  5.5b  is  zero,  the  solution  for 
f}f  is  readily  obtained  and  reads 


fn 


fb™  exp 


h 

eFn 


/  Y ~]Sop(k,k')\kx=z  +  i(qvF  ~  nui) 

J  k°  1  - 


dz 


(5.7) 


where  the  function  /blra  is  the  solution  fjy  evaluated  at  the  boundary  kx  =  k°  =  sj ( kc )2  —  ( ky )2.  On 
substituting  Eq.  5.7  into  Eq.  5.5a,  and  evaluating  the  resulting  equation  at  the  boundary  leads  to  a 
homogeneous  integral  equation  of  the  form, 

fbn(ky)=  r  dyfln(y)K(ky,  y,  n),  (5.8) 

J -k° 

where  K  is  a  complex  kernel.  Eq.  5.8  can  be  solved  numerically  by  discretizing  the  integral,76  which  leads  to 
a  unitary  matrix  equation  of  the  form  K  f  =  /.  It  is  found  that  these  eigenvectors  exist  only  for  qvF  =  nw, 
which  reduces  Eq.  5.5  to  the  steady  state-homogenous  equation  for  fh  =  f ft.64  A  similar  analysis  on  the 
harmonics  /ml,m  >  1  shows  that  solutions  exist  only,  for  mqvF  —  w,  so  fml  also  reduces  to  /q.  Since  fmn 
are  related  to  jO-LO-1)  by  Eq.  5.5,  it  can  be  deduced  from  Eq.  5.5  that  solutions  for  all  the  harmonics 
fmn  sucb  that  m  ^  n  >  0  that  exist  only  when  ui  =  mqvp/n ,  yield  fmn  =  /q  ,  and  should  therefore  be 
discarded.  As  a  result,  only  the  harmonics  fmn  such  that  m,n  >  0  and  m  =  n  >  0  are  the  only  remaining 
terms  of  the  series  (5.4)  for  f(k,x,t). 

The  current  density  on  the  plane  is  given  by, 


Jx{x ,  t)  =  —4  evF  E  x,  t )  cos  {(j)') 

k' 


(5.9) 


where  </>'  is  the  angle  between  k'  and  the  k'x  axis,  and  5ft/(fc',x,t)  is  the  real  part  of  the  distribution  func¬ 
tion.  In  this  analysis,  the  a-c  frequency  is  written  in  units  of  uj/lof,  where  u>f  =  2ttcF0vf /kioop,  and  the 
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Figure  5.2:  (Color  online)  (a)  Cross-sections  of  real  parts  of  the  harmonics  /n,/22,/33,  and  homogeneous  steady 
state  distribution  function  fh.  (b)Imaginary  parts  and  fh.  Dash  (40X/11),  dot  (100 Xf22),  dash-dot  (100 Xf33), 
and  solid  ( fh ).  Fa  =  1  kV/cm,  F\/F0  =  0.1,  uj/ujf  =  0.5,  td/ti  =  0.5,  and  7  =  0.01. 


wavelength  is  modulated  by  a  dimensionless  parameter  a  such  that 


q  = 


2i ra 
Xf 


(5.10) 


where  A f  =  kujop/eF0,  and  0  <  a  <  1.  As  in  a  previous  work,  '7  the  dimensionless  damping  parameter  7 
given  by 


7  = 


(5.11) 


To  Sop{k,  fc')U=1.5fec 

is  used  to  modulate  the  strength  of  low  energy  collisions  compared  to  OP  scattering. 

Figure  5.2  shows  cross  sections  of  the  real  parts  (5.2a),  and  imaginary  parts  (5.2b)  of  the  first  three 
harmonics  /n,/22,/33,  with  the  steady  state  homogeneous  distribution  fh  along  ky  =  0.  The  applied  field 
in  this  case  is  Fa  =  lkV/cm,  such  that  F\/F0  =  0.1,  w/wf  =  0.5,  and  a  =  1.  The  low  energy  damping 
parameter  is  set  at  7  =  0.01,  and  t0/ti  =  0.5.  From  the  figure,  it  can  be  seen  that  the  harmonics  oscillate  as 
a  function  of  kx  with  amplitude  fmm  decreasing  as  m  increases,  justifying  the  choice  of  the  solution  (5.4). 
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Also,  note  that  /  is  always  positive  as  fh  3>  fmn. 


Figure  5.3:  (Color  online)  Current  density  amplitude  versus  frequency  for  different  values  of  7  and  F0.  (a)  First 
harmonic,  (b)  Second  harmonic,  (c)  Third  harmonic.  Top  row  (7  =  0.01),  middle  row  (7  =  0.1),  and  bottom  row 
(7  =  0.3).  Solid  (F0  =  0.5  kV/cm),  dash  ( Fa  =  1  kV/cm),  and  dot  (F0  =  1.5  kV/cm).  F\/F0  =  0.1,  r0/n  =  0.5,  and 
a  —  1. 


Figure  5.3  displays  the  current  density  amplitude  relative  to  the  first  three  harmonics  as  a  function  of 
frequency  for  different  values  of  Fa  and  7.  As  in  Fig.  5.2,  the  applied  field  is  such  that  Fi/F0  =  0.1, 
t0/ti  =  0.5,  and  a  =  1.  For  all  three  harmonics,  there  are  resonances  that  occur  at  u  =  ljf,  for  which  the 
amplitudes  fmm(kx,ky)  are  real,  and  the  currents  are  in  phase  with  the  applied  field.  By  considering  the 
first  row  of  Fig.  5.3  (7  =  0.01),  one  observes  that  there  is  no  resonance  for  F0  =  0.5kV/cm  (first  column). 
This  is  because  the  contribution  of  the  negative  values  of  f11(kx,ky)  in  the  total  current  (Eq.  5.9)  offset 
that  of  the  positive  values  at  w  =  uj?.  As  the  field  increases  from  0.5kV/cm  to  lkV/cm,  more  electrons 
are  able  to  escape  low  energy  scattering  in  region  I  and  reach  the  boundary  k  =  kc  and  beyond,  so  the 
amplitude  fll(kx,  ky)  increases  in  region  II,  inducing  a  net  positive  current,  so  as  to  achieve  resonance.  For 
all  three  harmonics,  the  resonance  Quality  factor,  Q  =  lof/ Aojfwhm ,  where  A lofwhm  is  the  full  width 
of  the  current  density  profile  at  half  maximum  value,  first  increases  ( Q  ss  24.9  to  Q  «  25.2  for  the  second, 
and  Q  m  23  to  Q  ss  28  for  the  third)  as  the  field  increases  from  0.5kV/cm  to  lkV/cm.  A  further  increase 
in  the  field  to  1.5kV/cm  causes  more  electrons  to  reach  region  II  and  scatter  with  OP’s,  broadening  the 
distribution  in  the  process.  As  a  result,  Q  decreases  (from  Q  ss  20  to  Q  ~  13  for  the  first  harmonic,  from 
Q  ss  25.2  to  Q  «  24  for  the  second,  and  from  Q  «  28  to  Q  «  27  for  the  third)  as  the  field  increase  from 
lkV/cm  to  1.5kV/cm. 
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As  low  energy  scattering  (7)  increases,  the  system  achieves  resonance  even  at  low  fields  (Fig.  5.3  first 
column).  For  the  first  harmonic,  an  increase  in  7  makes  the  amplitude  more  positive  in  region  J,  and  the 
overall  current  becomes  positive,  achieving  resonance.  One  can  also  see  that  the  current  density  amplitude 
increases,  while  the  overall  resonance  Q-  factor  decreases  as  low  energy  scattering  increases.  For  the  first 
harmonic  with  7  =  0.1  (second  row,  first  column),  one  get  Q  «  0.7  (Fa  =  0.5kV/cm),  Q  ks  0.7  (F0  = 
lkV/cm),  and  6  ( F0  =  1.5kV/cm),  only.  This  is  due  to  the  fact  that  low  energy  scattering  redistribute 
carriers  towards  the  k  =  0  region.  As  seen  in  a  previous  work,68  carrier  interactions  with  OP’s  are  essential 
for  resonance  to  occur.  Higher  7-damping  causes  fewer  carriers  to  interact  with  OP’s  and  the  resonance 
Q-factor  decreases.  The  current  density  amplitude  increases  because  the  distribution  increases  around  k  =  0. 
Note  that  when  ui  =  u;?,  the  system  is  a  mixture  of  modes  of  oscillations  with  resonant  frequencies  mujF- 


J 


22 


CO/COp 


Figure  5.4:  (Color  online)  Second  harmonic  current  density  amplitude  vs  frequency  for  different  values  of  a.  Solid 
(a  =  0.8),  dash  (a  =  1.0)  and  dot  (a  =  1.5).  F\/F0  =  0.1,  t0/ti  =  0.5,  and  7  =  0.01. 
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Figure  5.4  shows  plots  of  the  second  harmonic  current  density  versus  frequency  for  different  values  of  a 
(gate  spacing).  The  damping  in  this  case  is  again  set  at  7  =  0.01,  F\/F0  =  0.1,  and  t0/t\  =  0.5.  From 
the  figure,  the  observed  resonances  are  achieved  at  the  frequencies  u>f,  0.8 wp,  and  1.5 u>f  corresponding  to 
a  =  1,  a  =  0.8,  and  a  =  1.5  respectively,  which  shows  that  for  a  particular  F„.  the  resonance  frequency 
is  tunable  with  the  wavelength  of  the  applied  field.  This  important  result  indicates  the  potential  use  of 
graphene  in  terahertz  sources  and  detectors.68,77 

In  conclusion,  an  analysis  of  carrier  dynamics  in  graphene  subjected  to  periodic,  time,  and  space  depen¬ 
dent  electric  fields  and  scattering  times  has  been  carried  out  in  this  chapter.  The  model  shows  that,  high  Q 
resonances  can  be  achieved  when  lo  =  u>f-  As  expected,  the  Q-factors  decrease  with  damping  7.  Another 
observation  is  that  at  resonance,  the  system  consists  of  carrier  excitations  with  frequencies  mojF ,  for  m  >  1. 
As  a  result  the  system  is  essentially  a  mixer  since  an  input  frequency  ljf,  creates  the  harmonics  mwp , 
and  appropriate  filters  should  be  used  to  pick  out  the  required  frequencies.  Also,  the  resonant  frequency  is 
tunable  with  the  wavelength  of  the  applied  field.  Note  that  the  wavelength  A  of  the  a-c  field  is  twice  the 
distance  between  successive  gates  in  Fig.  5.1.  Consequently,  graphene  can  potentially  be  used  to  make  high 
power,  tunable  terahertz  devices  that  operate  at  room  temperature. 
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Chapter  6 


CONCLUSIONS 


This  dissertation  presents  a  study  of  the  dynamics  of  carriers  in  graphene  under  the  influence  of  external 
electric  fields.  In  particular,  this  dissertation  analyses  the  interplay  between  the  electric  fields  and  OP 
scattering  to  generate  current  oscillations  with  frequencies  in  the  THz  range.  Such  oscillations  had  been 
predicted  and  observed  in  conventional  semiconductors.41,42  Room  temperature  observation  of  current  os¬ 
cillations  in  these  materials  is  however  not  possible.  This  is  because  at  room  temperature,  the  OP  energy 
( huiop  ~  0.04  eV)45  is  comparable  to  thermal  broadening  of  the  charge  distribution  and  the  oscillations  are 
immediately  damped.41,43,44  At  low  temperatures,  impurity  scattering  becomes  dominant  and  the  oscilla¬ 
tions  are  strongly  damped.  Impurity  scattering  can  be  reduced  by  lowering  the  dopant  density  which  in  turn 
lowers  the  carrier  density  reducing  the  oscillation  amplitude  in  the  process.  The  high  values  of  conductivity 
and  OP  energy  (hu>op  ~  0.2  eV) 38,39  of  graphene  make  it  possible  for  observing  room  temperature  THz 
oscillations. 

The  onset  of  current  oscillations  in  graphene  is  however  restricted  by  low  energy  scattering.  In  device 
applications,  charge  impurity  scattering  can  reduced  by  via  screening.  This  cab  be  achieved  by  placing  a 
graphene  layer  between  two  high  K  dielectrics.57  Screening  might  however  be  offset  by  interface  and  remote 
static  charges  contained  in  the  dielectrics. 58  Also,  the  presence  of  the  substrate  causes  carriers  in  graphene  to 
interact  with  remote  interface  phonons  (RIP),  increasing  low  energy  scattering.48,59,60,61  Freely  suspended 
graphene  layers  should  therefore  be  used  for  best  results. 

For  carrier  oscillations  to  occur  in  graphene,  the  applied  field  must  be  high  enough  that  carriers  escape 
low  energy  scattering  and  gain  enough  energy  to  interact  with  OP’s.  As  a  result,  the  process  requires  external 
biases  V  such  that  V  >  Hcoop/e.  The  field  should  not  be  too  high  so  that  carriers  are  almost  immediately 
scattered  by  OP’s  once  they  reach  energies  E  >  huiop ,  so  as  to  maintain  coherence  of  the  carrier  distribution 
function.  This  condition  is  satisfied  for  Tgr  top,  which  imposes  an  upper  limit  on  the  applied  field,  F  <C  5 
kV/cm  for  top  ~  0.4  ps.  These  conditions  also  require  that  sample  lengths  L  such  that  L  >  vftop ■  Because 
the  oscillations  are  damped  after  a  few  cycles,  the  sample  length  should  not  be  higher  than  a  few  values  of 
A  =  VpTgr. 
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This  dissertation  also  examines  the  effects  of  an  additional  spatially  uniform  ac  field  on  the  current 
oscillations  in  graphene.  As  discussed  earlier,  the  frequency  of  the  current  oscillations  depends  on  the 
applied  dc  field.  The  frequency  of  the  ac  component  modulates  the  frequency  of  carrier  oscillations,  and 
the  carriers  behave  as  a  parametric  oscillator.  Because  of  strongly  dissipative  nature  of  OP  scattering, 
the  system  considered  in  this  dissertation  behaves  differently  to  the  usual  parametric  oscillators.  Normal 
parametric  resonance  is  achieved  when  the  input  frequency  is  twice  the  frequency  of  natural  oscillations.67 
For  the  system  of  electrons  in  graphene,  resonance  is  achieved  for  input  frequencies  u  such  that  w  =  0.56wf, 
where  up  =  2n/Tgr.  When  only  the  dc  field  is  applied,  the  period  of  current  oscillations  is  t*  «  1.8 Tgr.  This 
is  because,  some  electrons  reach  energies  E  >  Hloop  and  are  not  immediately  scattered  by  OP’s.  These 
electrons  penetrate  deeper  into  the  high  energy  region,  stretching  the  period  of  the  oscillations.  As  a  result, 
the  frequency  of  such  oscillations  is  0.56wf,  thus  resonance  occurs  when  the  input  frequency  matches  the 
frequency  of  oscillations  with  just  the  dc  field  applied.  The  softness  and  anomaly  of  this  type  of  parametric 
resonance  can  be  attributed  to  the  effectiveness  of  OP  scattering  of  carriers  in  graphene.  The  same  analysis 
can  be  applied  to  a  system  of  electrons  under  the  influence  of  a  spatially  dependent  field  that  is  periodically 
modulated,  F( x)  =  F(x  +  d ),  under  steady  state  conditions, in  the  streaming  case.  The  resonance  condition 
in  that  case,  Fa  =  Hu;op/0.56ed,  can  be  the  basis  of  a  field  detector. 

Based  on  the  results  of  this  dissertation,  graphene  may  be  used  in  THz  devices  that  operate  at  room 
temperature.  The  imposed  lower  and  upper  limit  on  sample  lengths  ensure  that  such  devices  will  be  compact. 
The  insensitivity  of  the  resonance  condition  u>  =  t/ujf,  r/  ss  0.56,  to  the  applied  dc  field  Fa  ensures  that  tunable 
resonant  frequencies  covering  most  of  the  THz  range  may  be  obtained.  Sharper  resonances  at  w  =  wj?  may 
be  obtained  in  a  graphene  based  FET  with  periodic  gating.  Harmonics  with  frequencies  towf,  to  >  1 
are  also  generated  in  this  case.  The  fundamental  frequency  ujf  is  also  tunable  with  the  gate  separation. 
Consequently,  high  powered  devices  that  create  and  detect  electromagnetic  radiation  spanning  the  THz  gap 
may  be  realized. 
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Appendix 

Derivation  of  Equation  3.4 


Boltzmann  Transport  Equation(BTE)  in  the  two  regions  reads: 


dfL(k,t )  eF(t)  dfL(k,t ) 

dt  h  dkx 


h(k ,  t)  -  f0{k) 
tle 


Y,S0P$tifH$,t) 

k' 


(A. la) 


dfnjk,  t)  eF(t)  dfH(k,t) 

dt  H  dkx 


-fH(k,t)J2Sop(k,k') 

k' 


By  using  the  substitution  n  =  kx  +  (3(t),  and  choosing 


P(t) 


F(s)ds , 


(A.lb) 


(■2) 


Eq.  lb  reduces  to 

\F{p- \k  -  kx))jjr  =  ~9H  E  S0p{k,  k')  (.3) 

where  gn{kx ,  kv\  k)  =  fn(kx,  ky ,  j3~l  {k  —  kx)) ,  and  /3_1  is  the  inverse  function  of  /?  so  that  /3_1/3(f)  =  t.  The 
general  solution  of  Eq.  3  above  is  of  the  form: 


9H{kx,kv;K )  =gH(k°,ky ; 


dp  top  1(p,ky)\ 
F[/3-i(K-p)]  /’ 


where  kx  =  \  jkc 


-  ky 


is  the  kx  value  at  the  boundary  k  =  kc,  and 


(•4) 


top  l{kxiky)  -  Y  SOP(k ,  k'). 


Using  the  substitution  (2)  above,  the  time  dependent  general  solution  of  Eq.  lb  is  then 


fH(kx,ky,t)  =  fb(kv,/3  {/3(t)  +  kx-  kx))M(kx,ky,t), 


(.5) 
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where  fb(ky,t)  =  fH{kx,ky,t)  is  the  distribution  function  evaluated  at  the  boundary,  and 


M(kx,ky,t)  =  exp{  -  -  / 
^  e  J  k 


dp  tqp  1(p,kv ) 

ko  ^[/3_1(/3(t)  +  kx  -p)\ 


}• 


The  substitution  (2)  reduces  Eq.  la  to; 


dgp  _ _ h _ r  _  gL  ~  fc 

dkx  l  tle 


(  _  +  y-  Sop{%^  k)fH{k',  p~\K  -  h 

f  T  T  R  ^ 


(.6) 


Eq.  6  is  a  first  order  partial  differential  equation  with  a  solution  of  the  form: 

h  rk*  {fo(p,ky)  +  J2k'SOP{k',k)  fH{k',h)}ex  P{^} 

gL(kx,  ky- K)  =  -  exp  {  -  P  [K  kx)}  /  dp - 


7~LE  J  —  fcO 


Fit  i) 


(•7) 


where  t\  =  (3  1(k  —  p ).  Using  Eq.  2,  the  time  dependent  solution  of  Eq.  la  is  given  by 


fL{kx,ky,t)  =  -ex  p{ - —  }  [  dp{/o(p,fcy) +VS'op(fc',fc)  /»(£',  t2)\ exp  {  —  }/F(t2),  (.8) 

e  Tpp  J_k  o  t  ^  kx=p  >  tle 


where  t%  =  P  l{P{t)  +  kx  —p).  On  evaluating  Eq.  8  at  the  boundary  kx  =  kxl  using  the  boundary  condition 
fb{ky ,  t)  =  ky,t)  =  ky,t),  and  Eq.  5  leads  to  an  integral  equation  of  the  form: 


fb{ky,t)  =  fb(ky,t)  + 


—yy  [  dp  I  dk'SOP{k':k)  fb(kv',  t")M{k' ,  t')  exp{- - -}/-FV')>  (.9) 

27T  (3  J — /-0  J  kx — p  7~LE 


where  /3(t')  =  Pit)  +  kx  —  p,  Pit”)  =  Pit')  +  kx  —  \J  1  —  ( ky ')2,  A:;c/  =  k’  cos(c/)'),  ky  =  fc'sin(0'),  and 


fb(ky,t)  =  -  I  dpf0{p,ky)ex p{- - 

e  y_fco  Tip 


(.10) 
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